1 Eukaryotic virome analysis

#### load phyloseq project ####
GP.M1 <- readRDS("~/Desktop/mosq_Guad_infection/RDS/GIM_virome_raw.RDS")
GP.M1
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 152 taxa and 154 samples ]
## sample_data() Sample Data:       [ 154 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 152 taxa by 7 taxonomic ranks ]

1.1 Eukaryotic virome of Aedes aegypti

#### color ####
library(viridis)
col = plasma(9,alpha = 1, begin=0.9, end = 0, direction = 1)
scales::show_col(col)

1.1.1 Subset eukaryotic virome of Aedes aegypti

GP.M1_AA <- subset_samples(GP.M1, mosq_species  == "Ae_aegypti")
GP.M1_AA <- subset_samples(GP.M1_AA, group.plaque.qPCR.1000. != "")
GP.M1_AA <- subset_taxa(GP.M1_AA, species != "Chuvirus Mos8Chu0")
GP.M1_AA <- subset_taxa(GP.M1_AA, species != "Kaiowa virus")
GP.M1_AA <- subset_taxa(GP.M1_AA, species != "Cumbaru virus")
GP.M1_AA <- subset_taxa(GP.M1_AA, species != "Guato virus")
GP.M1_AA <- subset_taxa(GP.M1_AA, species != "Trichoplusia ni TED virus")
GP.M1_AA <- subset_taxa(GP.M1_AA, species != "Guadeloupe Culex tymo-like virus")
GP.M1_AA
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 126 taxa and 92 samples ]
## sample_data() Sample Data:       [ 92 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 126 taxa by 7 taxonomic ranks ]
sample_data(GP.M1_AA)$Treatment <- factor(sample_data(GP.M1_AA)$Treatment, levels = c("AA-sucrose", "AA-blood", "AA-ZIKV"))
sample_data(GP.M1_AA)$Info1 <- factor(sample_data(GP.M1_AA)$Info1, levels = c("AA-7dpi-sucrose", "AA-7dpi-blood", "AA-7dpi-ZIKV", "AA-21dpi-sucrose","AA-21dpi-blood","AA-21dpi-ZIKV"))
sample_data(GP.M1_AA)$group.plaque.qPCR.0.body.0. <- factor(sample_data(GP.M1_AA)$group.plaque.qPCR.0.body.0, levels = c("AA-7dpi-sucrose", "AA-7dpi-blood", "AA-7dpi-ZIKV(-/-/-)", "AA-7dpi-ZIKV(-/-/+)", "AA-7dpi-ZIKV(-/+/+)", "AA-21dpi-sucrose", "AA-21dpi-blood", "AA-21dpi-ZIKV(-/-/-)", "AA-21dpi-ZIKV(-/-/+)",  "AA-21dpi-ZIKV(-/+/-)", "AA-21dpi-ZIKV(-/+/+)", "AA-21dpi-ZIKV(+/+/+)"))
sample_data(GP.M1_AA)$plaque_assay <- factor(sample_data(GP.M1_AA)$plaque_assay,
                                             levels = c("AA-7dpi-sucrose", "AA-7dpi-blood", "AA-7dpi-ZIKV(-)",
                                                        "AA-21dpi-sucrose", "AA-21dpi-blood", "AA-21dpi-ZIKV(-)", "AA-21dpi-ZIKV(+)"))

sample_data(GP.M1_AA)$dpi <- factor(sample_data(GP.M1_AA)$dpi, levels = c("7dpi", "21dpi"))

1.1.2 Merge on species level

GP.M1_species <- GP.M1_AA %>% tax_glom(taxrank = "species")
GP.M1_species
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 55 taxa and 92 samples ]
## sample_data() Sample Data:       [ 92 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 55 taxa by 7 taxonomic ranks ]
GP.M1_samselect = prune_samples(sample_sums(GP.M1_species) >0, GP.M1_species)

# species
GP.M1_species_select_0 = prune_taxa(taxa_sums(GP.M1_samselect) > 0, GP.M1_samselect)
GP.M1_species_select_0
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 40 taxa and 92 samples ]
## sample_data() Sample Data:       [ 92 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 40 taxa by 7 taxonomic ranks ]
GP.M1_species_select = prune_taxa(taxa_sums(GP.M1_samselect) > 500, GP.M1_samselect)
GP.M1_species_select
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 19 taxa and 92 samples ]
## sample_data() Sample Data:       [ 92 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 19 taxa by 7 taxonomic ranks ]

1.1.3 Heatmap on viral species level

#### get table of reads for each viral species ####
OTU_species <- merge(otu_table(GP.M1_species_select), tax_table(GP.M1_species_select), by=0)
rownames(OTU_species) <- OTU_species$species

#### order viral species by familiy #### 
OTU_species$family <- as.character(OTU_species$family)
OTU_species$family
##  [1] "Virgaviridae"     "Rhabdoviridae"    "Flaviviridae"     "Orthomyxoviridae"
##  [5] "Flaviviridae"     "Totiviridae"      NA                 "Xinmoviridae"    
##  [9] "Totiviridae"      NA                 "Totiviridae"      "Flaviviridae"    
## [13] NA                 "Reoviridae"       "Flaviviridae"     NA                
## [17] "Phenuiviridae"    "Flaviviridae"     "Xinmoviridae"
OTU_species$family[is.na(OTU_species$family)] <- "unclassified"
OTU_species_m<- melt(OTU_species)
## Using Row.names, superkingdom, phylum, class, order, family, genus, species as id variables
OTU_species_m_agg <- aggregate(value~family+species, FUN=sum, data=OTU_species_m)
OTU_species_m_agg1 <- aggregate(value~family, FUN=sum, data=OTU_species_m)
family_ord <- OTU_species_m_agg1[order(-OTU_species_m_agg1$value),]$family
family_ord <- c(family_ord[-3],"unclassified")
family_ord
## [1] "Phenuiviridae"    "Totiviridae"      "Flaviviridae"     "Orthomyxoviridae"
## [5] "Xinmoviridae"     "Virgaviridae"     "Reoviridae"       "Rhabdoviridae"   
## [9] "unclassified"
OTU_species_m_agg$index <- 1
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Totiviridae")] <- 2
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Flaviviridae")] <- 3
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Orthomyxoviridae")] <- 4
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Xinmoviridae")] <- 5
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Metaviridae")] <- 6
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Virgaviridae")] <- 7
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Rhabdoviridae")] <- 8
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Reoviridae")] <- 9
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "unclassified")] <- 10
OTU_species_m_agg <- OTU_species_m_agg[order(OTU_species_m_agg$index, -OTU_species_m_agg$value),]
OTU_species_m_agg
##              family                                   species   value index
## 14    Phenuiviridae             Phasi Charoen-like phasivirus 8253667     1
## 2       Totiviridae                   Aedes aegypti totivirus 7849526     2
## 1       Totiviridae             Aedes aegypti toti-like virus   43484     2
## 5       Totiviridae            Australian Anopheles totivirus    1722     2
## 19     Flaviviridae                                Zika virus  409642     3
## 8      Flaviviridae                   Cell fusing agent virus   55036     3
## 16     Flaviviridae                 unclassified Flaviviridae    8111     3
## 13     Flaviviridae                        Menghai flavivirus    6987     3
## 18     Flaviviridae            Xishuangbanna aedes flavivirus     915     3
## 9  Orthomyxoviridae Guadeloupe mosquito quaranja-like virus 1   57556     4
## 4      Xinmoviridae                          Aedes anphevirus   21985     5
## 15     Xinmoviridae                   unclassified Anphevirus   18618     5
## 7      Virgaviridae              Bombus-associated virus Vir1    1644     7
## 17    Rhabdoviridae                unclassified Rhabdoviridae     741     8
## 12       Reoviridae                           Liao ning virus     755     9
## 10     unclassified                 Guadeloupe mosquito virus  583366    10
## 6      unclassified                  Beihai barnacle virus 12  231777    10
## 11     unclassified                  Humaita-Tubiacanga virus    7435    10
## 3      unclassified  Aedes alboannulatus orthomyxo-like virus     868    10
species_ord <- as.character(OTU_species_m_agg$species)
species_ord
##  [1] "Phasi Charoen-like phasivirus"            
##  [2] "Aedes aegypti totivirus"                  
##  [3] "Aedes aegypti toti-like virus"            
##  [4] "Australian Anopheles totivirus"           
##  [5] "Zika virus"                               
##  [6] "Cell fusing agent virus"                  
##  [7] "unclassified Flaviviridae"                
##  [8] "Menghai flavivirus"                       
##  [9] "Xishuangbanna aedes flavivirus"           
## [10] "Guadeloupe mosquito quaranja-like virus 1"
## [11] "Aedes anphevirus"                         
## [12] "unclassified Anphevirus"                  
## [13] "Bombus-associated virus Vir1"             
## [14] "unclassified Rhabdoviridae"               
## [15] "Liao ning virus"                          
## [16] "Guadeloupe mosquito virus"                
## [17] "Beihai barnacle virus 12"                 
## [18] "Humaita-Tubiacanga virus"                 
## [19] "Aedes alboannulatus orthomyxo-like virus"
which(species_ord == "Zika virus")
## [1] 5
species_ord_1 <- c("Zika virus", species_ord[-5])

#### transform the table - reads for each viral species ##### 
rm <- c("Row.names", "superkingdom", "phylum", "class", "order", "family", "genus", "species")
OTU_species <- OTU_species[,!colnames(OTU_species) %in% rm]
OTU_species_trans <- log10(OTU_species)
is.na(OTU_species_trans) <- sapply(OTU_species_trans, is.infinite)
OTU_species_trans[is.na(OTU_species_trans)]<-0
df.original <- OTU_species
df <- OTU_species_trans

#### use meta data for top annotation ####
meta = read.csv("~/Desktop/mosq_Guad_infection/GIM_metadata_new.csv", header=TRUE, row.names=1, sep=",")
meta_s = meta[!meta$mosq_species == "",]

meta_AA <- meta[colnames(df),]

#### order sample by ZIKV reads within each group & column annotation  #### 
meta_AAz <- merge(meta_AA, as.data.frame(t(df)),by=0)
head(meta_AAz[1:5, 1:5])
##   Row.names         mosq    body      mosq_or   head
## 1    GIM100 AaZ-21dpi-20 GIM-100 AMT-21dpi-20 GIM-80
## 2     GIM22   AaZ-7dpi-2  GIM-22   AMT-7dpi-2  GIM-2
## 3     GIM23   AaZ-7dpi-3  GIM-23   AMT-7dpi-3  GIM-3
## 4     GIM24   AaZ-7dpi-4  GIM-24   AMT-7dpi-4  GIM-4
## 5     GIM25   AaZ-7dpi-5  GIM-25   AMT-7dpi-5  GIM-5
meta_AAz$Treatment <- as.character(meta_AAz$Treatment)
meta_AAz$Treatment <- factor(meta_AAz$Treatment, levels = c("AA-sucrose", "AA-blood", "AA-ZIKV"))
meta_AAz$Info1 <- as.character(meta_AAz$Info1)
meta_AAz$Info1 <- factor(meta_AAz$Info1, levels = c("AA-7dpi-sucrose", "AA-7dpi-blood", "AA-7dpi-ZIKV", "AA-21dpi-sucrose","AA-21dpi-blood","AA-21dpi-ZIKV"))

meta_AAz$qPCR_head._ZIKV_WNV <- log10(meta_AAz$qPCR_head._ZIKV_WNV)
meta_AAz$qPCR_body._ZIKV_WNV <- log10(meta_AAz$qPCR_body._ZIKV_WNV)

which(colnames(meta_AAz)=="qPCR_head._ZIKV_WNV")
## [1] 16
which(colnames(meta_AAz)=="qPCR_body._ZIKV_WNV")
## [1] 17
is.na(meta_AAz[,16:17])<-sapply(meta_AAz[,16:17], is.infinite)
meta_AAz[,16:17][is.na(meta_AAz[,16:17])]<-0

meta_AAz$group_n1 <- meta_AAz$plaque_assay
meta_AAz$plaque_assay <- factor(meta_AAz$plaque_assay, levels = c("AA-7dpi-sucrose", "AA-7dpi-blood", "AA-7dpi-ZIKV(-)", "AA-21dpi-sucrose", "AA-21dpi-blood", "AA-21dpi-ZIKV(-)","AA-21dpi-ZIKV(+)"), ordered = TRUE)

meta_AAz_order1 <- meta_AAz[meta_AAz$dpi== "7dpi",][with(meta_AAz[meta_AAz$dpi== "7dpi",], order(plaque_assay, -qPCR_head._ZIKV_WNV, -qPCR_body._ZIKV_WNV, -`Zika virus`, -`Guadeloupe mosquito virus`,-`Phasi Charoen-like phasivirus`, -`Aedes aegypti toti-like virus`,-`Aedes anphevirus`,-`Guadeloupe mosquito quaranja-like virus 1`)), ]
meta_AAz_order2 <- meta_AAz[meta_AAz$dpi== "21dpi",][with(meta_AAz[meta_AAz$dpi== "21dpi",], order(plaque_assay, -qPCR_head._ZIKV_WNV,  -qPCR_body._ZIKV_WNV, -`Zika virus`, -`Guadeloupe mosquito virus`,-`Phasi Charoen-like phasivirus`,-`Aedes aegypti toti-like virus`,-`Aedes anphevirus`,-`Guadeloupe mosquito quaranja-like virus 1`)), ]
meta_AAz_order <- rbind(meta_AAz_order1,meta_AAz_order2)

df <- df[,meta_AAz_order$Row.names]
col_ha = HeatmapAnnotation(qPCR = anno_lines(meta_AAz_order[,16:17], gp = gpar(col = c("black","darkgrey")), height = unit(1.5, "cm"),ylim = c(-0.5,8), axis_param = list(at = c(0,2,4,6,8)),add_points = TRUE, pt_gp = gpar(col = c("black","darkgrey")), pch = c(16, 17)))

#### column split ####
df_1 <- df
df_1$ZIKV <- "others"
df_1$ZIKV[which(rownames(df_1)=="Zika virus")] <- "aZIKV"
#### draw figure #####
heatmap_col_tp2 <- colorRampPalette(brewer.pal(9, "YlGnBu"), space = "rgb")(100)
Heatmap(as.matrix(df),
        name = "log10 reads number", #title of legend
        column_title = "Samples", row_title = "viral species",
        row_names_gp = gpar(fontsize = 10),
        column_names_gp = gpar(fontsize = 6),
        column_title_gp = gpar(fontsize = 12),
        cluster_rows = F,
        cluster_columns= F,
        row_order = species_ord_1,
        column_labels = meta_AAz_order$mosq,
        top_annotation = col_ha,
        col = heatmap_col_tp2,
        column_split = meta_AAz_order$group_n1,
        row_split = df_1$ZIKV,
        border = TRUE,
        rect_gp = gpar(col = "white", lwd = 0.5))

1.1.4 Alpha diversity on species level

plot_richness <- function (physeq, x = "samples", color = NULL, shape = NULL, 
                           title = NULL, scales = "free_y", nrow = 1, shsi = NULL, measures = NULL, 
                           sortby = NULL) 
{
  erDF = estimate_richness(physeq, split = TRUE, measures = measures)
  measures = colnames(erDF)
  ses = colnames(erDF)[grep("^se\\.", colnames(erDF))]
  measures = measures[!measures %in% ses]
  if (!is.null(sample_data(physeq, errorIfNULL = FALSE))) {
    DF <- data.frame(erDF, sample_data(physeq))
  }
  else {
    DF <- data.frame(erDF)
  }
  if (!"samples" %in% colnames(DF)) {
    DF$samples <- sample_names(physeq)
  }
  if (!is.null(x)) {
    if (x %in% c("sample", "samples", "sample_names", "sample.names")) {
      x <- "samples"
    }
  }
  else {
    x <- "samples"
  }
  mdf = reshape2::melt(DF, measure.vars = measures)
  mdf$se <- NA_integer_
  if (length(ses) > 0) {
    selabs = ses
    names(selabs) <- substr(selabs, 4, 100)
    substr(names(selabs), 1, 1) <- toupper(substr(names(selabs), 
                                                  1, 1))
    mdf$wse <- sapply(as.character(mdf$variable), function(i, 
                                                           selabs) {
      selabs[i]
    }, selabs)
    for (i in 1:nrow(mdf)) {
      if (!is.na(mdf[i, "wse"])) {
        mdf[i, "se"] <- mdf[i, (mdf[i, "wse"])]
      }
    }
    mdf <- mdf[, -which(colnames(mdf) %in% c(selabs, "wse"))]
  }
  if (!is.null(measures)) {
    if (any(measures %in% as.character(mdf$variable))) {
      mdf <- mdf[as.character(mdf$variable) %in% measures, 
      ]
    }
    else {
      warning("Argument to `measures` not supported. All alpha-diversity measures (should be) included in plot.")
    }
  }
  if (!is.null(shsi)) {
    warning("shsi no longer supported option in plot_richness. Please use `measures` instead")
  }
  if (!is.null(sortby)) {
    if (!all(sortby %in% levels(mdf$variable))) {
      warning("`sortby` argument not among `measures`. Ignored.")
    }
    if (!is.discrete(mdf[, x])) {
      warning("`sortby` argument provided, but `x` not a discrete variable. `sortby` is ignored.")
    }
    if (all(sortby %in% levels(mdf$variable)) & is.discrete(mdf[, 
                                                                x])) {
      wh.sortby = which(mdf$variable %in% sortby)
      mdf[, x] <- factor(mdf[, x], levels = names(sort(tapply(X = mdf[wh.sortby, 
                                                                      "value"], INDEX = mdf[wh.sortby, x], mean, na.rm = TRUE, 
                                                              simplify = TRUE))))
    }
  }
  richness_map = aes_string(x = x, y = "value", colour = color, 
                            shape = shape)
  p = ggplot(mdf, richness_map) + geom_boxplot(na.rm = TRUE,outlier.shape = NA)+
    # geom_point(size=1)
    geom_jitter(position=position_jitter(0.25), size = 3.5, alpha=0.7)
  # if (any(!is.na(mdf[, "se"]))) {
  #   p = p + geom_errorbar(aes(ymax = value + se, ymin = value - 
  #                               se), width = 0.1)
  # }
  # p = p + theme(axis.text.x = element_text(angle = -90, vjust = 0.5, 
                                           # hjust = 0))
  p = p + ylab("Alpha Diversity Measure")
  p = p + facet_wrap(~variable, nrow = nrow, scales = scales)
  if (!is.null(title)) {
    p <- p + ggtitle(title)
  }
  return(p)
}

#### on diet ####
plot_richness(GP.M1_species_select_0, x="Info1",measures=c("Shannon"),col="Treatment", shape = "dpi") + 
  scale_color_manual(values = rev(c(col[9],col[5],col[2]))) +
  scale_shape_manual(values = c(16, 17))+
  ylim(-0.001, 2.5)+
  theme_bw()+
  theme(axis.title.x = element_blank(),
        axis.text.x  = element_text(size = 12, colour = "black", angle = 90, hjust = 1),
        strip.text.x = element_text(size = 12, colour = "black"),
        plot.title = element_text(hjust = 0.5, size = 12)) +
  ggtitle("Alpha diversity of eukaryotic virus in Ae. aegypti")+
  ylab("Alpha Diversity Value")

erich <- estimate_richness(GP.M1_species_select_0, measures = c("Shannon"))
pairwise.wilcox.test(erich$Shannon,sample_data(GP.M1_species_select_0)$Info1, p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum exact test 
## 
## data:  erich$Shannon and sample_data(GP.M1_species_select_0)$Info1 
## 
##                  AA-7dpi-sucrose AA-7dpi-blood AA-7dpi-ZIKV AA-21dpi-sucrose
## AA-7dpi-blood    0.57359         -             -            -               
## AA-7dpi-ZIKV     4.1e-05         0.03450       -            -               
## AA-21dpi-sucrose 0.02381         1.00000       0.00437      -               
## AA-21dpi-blood   0.03450         0.68452       0.24051      0.73980         
## AA-21dpi-ZIKV    0.00092         0.57359       4.1e-05      0.41492         
##                  AA-21dpi-blood
## AA-7dpi-blood    -             
## AA-7dpi-ZIKV     -             
## AA-21dpi-sucrose -             
## AA-21dpi-blood   -             
## AA-21dpi-ZIKV    0.70485       
## 
## P value adjustment method: BH
#### on diet and infection status ####
sample_data(GP.M1_species_select_0)$hd <- as.character(sample_data(GP.M1_species_select_0)$group.plaque.qPCR.1000.body.1000.)

sample_data(GP.M1_species_select_0)$hd[sample_data(GP.M1_species_select_0)$group.plaque.qPCR.1000.body.1000. == "AA-7dpi-ZIKV(-/+/+)"] <- "7dpe-Head(+)/Body(+)"
sample_data(GP.M1_species_select_0)$hd[sample_data(GP.M1_species_select_0)$group.plaque.qPCR.1000.body.1000. == "AA-7dpi-ZIKV(-/-/+)"] <- "7dpe-Head(-)/Body(+)"
sample_data(GP.M1_species_select_0)$hd[sample_data(GP.M1_species_select_0)$group.plaque.qPCR.1000.body.1000. == "AA-7dpi-ZIKV(-/-/-)"] <- "7dpe-Head(-)/Body(-)"
sample_data(GP.M1_species_select_0)$hd[sample_data(GP.M1_species_select_0)$group.plaque.qPCR.1000.body.1000. == "AA-21dpi-ZIKV(+/+/+)"] <- "21dpe-Head(+)/Body(+)"
sample_data(GP.M1_species_select_0)$hd[sample_data(GP.M1_species_select_0)$group.plaque.qPCR.1000.body.1000. == "AA-21dpi-ZIKV(-/+/+)"] <- "21dpe-Head(+)/Body(+)"
sample_data(GP.M1_species_select_0)$hd[sample_data(GP.M1_species_select_0)$group.plaque.qPCR.1000.body.1000. == "AA-21dpi-ZIKV(-/-/+)"] <- "21dpe-Head(-)/Body(+)"
sample_data(GP.M1_species_select_0)$hd[sample_data(GP.M1_species_select_0)$group.plaque.qPCR.1000.body.1000. == "AA-21dpi-ZIKV(-/-/-)"] <- "21dpe-Head(-)/Body(-)"

sample_data(GP.M1_species_select_0)$hd <- factor(sample_data(GP.M1_species_select_0)$hd,
                                             levels = c("AA-7dpi-sucrose","AA-7dpi-blood","7dpe-Head(-)/Body(-)","7dpe-Head(-)/Body(+)", "7dpe-Head(+)/Body(+)", 
                                                        "AA-21dpi-sucrose","AA-21dpi-blood","21dpe-Head(-)/Body(-)","21dpe-Head(-)/Body(+)", "21dpe-Head(+)/Body(+)"))

plot_richness(GP.M1_species_select_0, x="hd", measures=c("Shannon"), col="hd", shape = "dpi") + 
  scale_color_manual(values = rev(c(col[9],col[7],col[5],col[3],col[1],col[9],col[7],col[5],col[3],col[1]))) +
  scale_shape_manual(values = c(16, 17))+
  ylim(-0.0001,2.5)+
  theme_bw()+
  theme(axis.title.x = element_blank(),
        axis.text.x  = element_text(size = 12, colour = "black", angle = 90, hjust = 1, vjust=0.5),
        strip.text.x = element_text(size = 12, colour = "black"),
        plot.title = element_text(hjust = 0.5, size = 12)) +
  ggtitle("Alpha diversity of eukaryotic virus in Ae. aegypti")+
  ylab("Alpha Diversity Value")

erich <- estimate_richness(GP.M1_species_select_0, measures = c("Shannon"))
pairwise.wilcox.test(erich$Shannon,sample_data(GP.M1_species_select_0)$hd, p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum exact test 
## 
## data:  erich$Shannon and sample_data(GP.M1_species_select_0)$hd 
## 
##                       AA-7dpi-sucrose AA-7dpi-blood 7dpe-Head(-)/Body(-)
## AA-7dpi-blood         0.6760          -             -                   
## 7dpe-Head(-)/Body(-)  0.0013          0.0887        -                   
## 7dpe-Head(-)/Body(+)  0.0112          0.1111        0.9481              
## 7dpe-Head(+)/Body(+)  0.0549          0.2500        0.7858              
## AA-21dpi-sucrose      0.0397          1.0000        0.0199              
## AA-21dpi-blood        0.0549          0.7766        0.2983              
## 21dpe-Head(-)/Body(-) 0.0199          0.6814        0.0149              
## 21dpe-Head(-)/Body(+) 0.5669          0.7895        1.0000              
## 21dpe-Head(+)/Body(+) 0.0013          0.6814        0.0054              
##                       7dpe-Head(-)/Body(+) 7dpe-Head(+)/Body(+)
## AA-7dpi-blood         -                    -                   
## 7dpe-Head(-)/Body(-)  -                    -                   
## 7dpe-Head(-)/Body(+)  -                    -                   
## 7dpe-Head(+)/Body(+)  0.7858               -                   
## AA-21dpi-sucrose      0.0511               0.0893              
## AA-21dpi-blood        0.3882               0.7766              
## 21dpe-Head(-)/Body(-) 0.0511               0.2036              
## 21dpe-Head(-)/Body(+) 1.0000               0.9000              
## 21dpe-Head(+)/Body(+) 0.0199               0.0624              
##                       AA-21dpi-sucrose AA-21dpi-blood 21dpe-Head(-)/Body(-)
## AA-7dpi-blood         -                -              -                    
## 7dpe-Head(-)/Body(-)  -                -              -                    
## 7dpe-Head(-)/Body(+)  -                -              -                    
## 7dpe-Head(+)/Body(+)  -                -              -                    
## AA-21dpi-sucrose      -                -              -                    
## AA-21dpi-blood        0.7967           -              -                    
## 21dpe-Head(-)/Body(-) 0.5669           0.7766         -                    
## 21dpe-Head(-)/Body(+) 0.5669           0.7895         0.5336               
## 21dpe-Head(+)/Body(+) 0.5436           0.7895         0.9160               
##                       21dpe-Head(-)/Body(+)
## AA-7dpi-blood         -                    
## 7dpe-Head(-)/Body(-)  -                    
## 7dpe-Head(-)/Body(+)  -                    
## 7dpe-Head(+)/Body(+)  -                    
## AA-21dpi-sucrose      -                    
## AA-21dpi-blood        -                    
## 21dpe-Head(-)/Body(-) -                    
## 21dpe-Head(-)/Body(+) -                    
## 21dpe-Head(+)/Body(+) 0.2493               
## 
## P value adjustment method: BH

1.1.5 PCoA on species level

#### 7 dpe ##### 
GP.M1_species_select_0_7dpi = subset_samples(GP.M1_species_select_0, dpi == "7dpi")
GP.M1_species_select_0_rmZIKV <- subset_taxa(GP.M1_species_select_0_7dpi, !taxa_names(GP.M1_species_select_0_7dpi) == "GIM628_NODE_1_length_10744_cov_806.867254")
GP.M1_species_select_0_rmZIKV
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 39 taxa and 42 samples ]
## sample_data() Sample Data:       [ 42 samples by 19 sample variables ]
## tax_table()   Taxonomy Table:    [ 39 taxa by 7 taxonomic ranks ]
sample_variables(GP.M1_species_select_0_rmZIKV)
##  [1] "mosq"                              "body"                             
##  [3] "mosq_or"                           "head"                             
##  [5] "Index"                             "mosq_species"                     
##  [7] "Treatment"                         "dpi"                              
##  [9] "group.plaque.qPCR.1000."           "group.plaque.qPCR.0."             
## [11] "group.plaque.qPCR.1000.body.1000." "group.plaque.qPCR.0.body.0."      
## [13] "Info1"                             "plaque_assay"                     
## [15] "qPCR_head._ZIKV_WNV"               "qPCR_body._ZIKV_WNV"              
## [17] "head_ZIKV.1000"                    "body_ZIKV.1000"                   
## [19] "hd"
GP.ord.GP.M1_species_select_0_rmZIKV <- ordinate(GP.M1_species_select_0_rmZIKV, "PCoA", "bray")

plot_ordination(GP.M1_species_select_0_rmZIKV, GP.ord.GP.M1_species_select_0_rmZIKV, type = "samples", color = "Treatment", shape = "group.plaque.qPCR.1000.body.1000.")+
  theme_bw()+
  geom_point(size = 4)  + 
  scale_shape_manual(values=rev(shap[2:6])) +
  scale_color_manual(values = c(col[2],col[5],col[9]))+
  theme(legend.title=element_blank(),
        legend.position = "right",
        plot.title = element_text(size = 12,hjust = 0.5))+
  ggtitle("PCoA of eukaryotic virus in Ae. aegypti at 7 dpe excluding ZIKV") +
  stat_ellipse(type = "norm", linetype = 2,level=0.6,aes(group = Treatment))

set.seed(1)
# Calculate bray curtis distance matrix
GP.M1_species_select_bray <- phyloseq::distance(GP.M1_species_select_0_rmZIKV, method = "bray")
# make a data frame from the sample_data
GP.M1_species_select_sampledf <- data.frame(sample_data(GP.M1_species_select_0_rmZIKV))
# Adonis test
adonis(GP.M1_species_select_bray ~ group.plaque.qPCR.1000.body.1000., data = GP.M1_species_select_sampledf)
## 
## Call:
## adonis(formula = GP.M1_species_select_bray ~ group.plaque.qPCR.1000.body.1000.,      data = GP.M1_species_select_sampledf) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                                   Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## group.plaque.qPCR.1000.body.1000.  4    3.2397 0.80992  2.5814 0.21818  0.001
## Residuals                         37   11.6087 0.31375         0.78182       
## Total                             41   14.8484                 1.00000       
##                                      
## group.plaque.qPCR.1000.body.1000. ***
## Residuals                            
## Total                                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#### 21 dpe ##### 
GP.M1_species_select_0_21dpi = subset_samples(GP.M1_species_select_0, dpi == "21dpi")
GP.M1_species_select_0_rmZIKV <- subset_taxa(GP.M1_species_select_0_21dpi, !taxa_names(GP.M1_species_select_0_21dpi) == "GIM628_NODE_1_length_10744_cov_806.867254")
GP.M1_species_select_0_rmZIKV
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 39 taxa and 50 samples ]
## sample_data() Sample Data:       [ 50 samples by 19 sample variables ]
## tax_table()   Taxonomy Table:    [ 39 taxa by 7 taxonomic ranks ]
GP.ord.GP.M1_species_select_0_rmZIKV <- ordinate(GP.M1_species_select_0_rmZIKV, "PCoA", "bray")

plot_ordination(GP.M1_species_select_0_rmZIKV, GP.ord.GP.M1_species_select_0_rmZIKV, type = "samples", color = "Treatment", shape = "group.plaque.qPCR.1000.body.1000.")+
  theme_bw()+
  geom_point(size = 4)  + 
  scale_shape_manual(values=c(shap[1:4],shap[6],shap[5])) +
  scale_color_manual(values = c(col[2],col[5],col[9]))+
  theme(legend.title=element_blank(),
        legend.position = "right",
        plot.title = element_text(size = 12,hjust = 0.5))+
  ggtitle("PCoA of eukaryotic virus in Ae. aegypti at 21 dpe excluding ZIKV")

set.seed(1)
# Calculate bray curtis distance matrix
GP.M1_species_select_bray <- phyloseq::distance(GP.M1_species_select_0_rmZIKV, method = "bray")
# make a data frame from the sample_data
GP.M1_species_select_sampledf <- data.frame(sample_data(GP.M1_species_select_0_rmZIKV))
# Adonis test
adonis(GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf)
## 
## Call:
## adonis(formula = GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## Treatment  2    0.7426 0.37129 0.96082 0.03928  0.473
## Residuals 47   18.1622 0.38643         0.96072       
## Total     49   18.9048                 1.00000

1.2 Eukaryotic virome of Culex quinquefasciatus

1.2.1 Subset eukaryotic virome of Culex quinquefasciatus

GP.M1_CQ <- subset_samples(GP.M1, mosq_species  == "Cx_quinquefasciatus")
GP.M1_CQ <- subset_taxa(GP.M1_CQ, species != "Kaiowa virus")
GP.M1_CQ
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 143 taxa and 60 samples ]
## sample_data() Sample Data:       [ 60 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 143 taxa by 7 taxonomic ranks ]

1.2.2 merge on species level

GP.M1_species <- GP.M1_CQ %>%
  tax_glom(taxrank = "species")
GP.M1_species
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 60 taxa and 60 samples ]
## sample_data() Sample Data:       [ 60 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 60 taxa by 7 taxonomic ranks ]
GP.M1_samselect = prune_samples(sample_sums(GP.M1_species) > 0, GP.M1_species)
GP.M1_samselect                                
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 60 taxa and 55 samples ]
## sample_data() Sample Data:       [ 55 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 60 taxa by 7 taxonomic ranks ]
GP.M1_species_select_0 = prune_taxa(taxa_sums(GP.M1_samselect) > 0, GP.M1_samselect)
GP.M1_species_select_0
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 27 taxa and 55 samples ]
## sample_data() Sample Data:       [ 55 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 27 taxa by 7 taxonomic ranks ]
GP.M1_species_select_0 <- subset_taxa(GP.M1_species_select_0, species != "Aedes aegypti totivirus")
GP.M1_species_select_0
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 26 taxa and 55 samples ]
## sample_data() Sample Data:       [ 55 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 26 taxa by 7 taxonomic ranks ]
GP.M1_species_select = prune_taxa(taxa_sums(GP.M1_samselect) > 500, GP.M1_samselect)
GP.M1_species_select <- subset_taxa(GP.M1_species_select, species != "Aedes aegypti totivirus")
GP.M1_species_select
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 9 taxa and 55 samples ]
## sample_data() Sample Data:       [ 55 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 9 taxa by 7 taxonomic ranks ]

1.2.3 Heatmap on viral species level

#### get table of reads for each viral species ####
OTU_species <- merge(otu_table(GP.M1_species_select), tax_table(GP.M1_species_select), by=0)
rownames(OTU_species) <- OTU_species$species
rownames(OTU_species)
## [1] "unclassified Chloriridovirus"      "West Nile virus"                  
## [3] "Lysoka partiti-like virus"         "unclassified Circoviridae"        
## [5] "Guadeloupe Culex tymo-like virus"  "Marine protobacilladnavirus 1"    
## [7] "Lampyris noctiluca errantivirus 1" "Wenzhou sobemo-like virus 3"      
## [9] "Phasi Charoen-like phasivirus"
#### order viral species by familiy #### 
OTU_species$family <- as.character(OTU_species$family)
OTU_species$family
## [1] "Iridoviridae"      "Flaviviridae"      "Picobirnaviridae" 
## [4] "Circoviridae"      NA                  "Bacilladnaviridae"
## [7] "Metaviridae"       NA                  "Phenuiviridae"
OTU_species$family[is.na(OTU_species$family)] <- "unclassified"
OTU_species_m<- melt(OTU_species)
## Using Row.names, superkingdom, phylum, class, order, family, genus, species as id variables
OTU_species_m_agg <- aggregate(value~family+species, FUN=sum, data=OTU_species_m)
OTU_species_m_agg1 <- aggregate(value~family, FUN=sum, data=OTU_species_m)
family_ord <- OTU_species_m_agg1[order(-OTU_species_m_agg1$value),]$family
family_ord <- c(family_ord[-1],"unclassified")
family_ord
## [1] "Flaviviridae"      "Phenuiviridae"     "Metaviridae"      
## [4] "Circoviridae"      "Picobirnaviridae"  "Iridoviridae"     
## [7] "Bacilladnaviridae" "unclassified"
OTU_species_m_agg$index <- 1
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Phenuiviridae")] <- 2
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Metaviridae")] <- 3
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Circoviridae")] <- 4
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Picobirnaviridae")] <- 5
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Iridoviridae")] <- 6
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Bacilladnaviridae")] <- 7
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "unclassified")] <- 8

OTU_species_m_agg <- OTU_species_m_agg[order(OTU_species_m_agg$index, -OTU_species_m_agg$value),]
OTU_species_m_agg
##              family                           species   value index
## 9      Flaviviridae                   West Nile virus   93203     1
## 5     Phenuiviridae     Phasi Charoen-like phasivirus    8886     2
## 2       Metaviridae Lampyris noctiluca errantivirus 1    1874     3
## 7      Circoviridae         unclassified Circoviridae    1597     4
## 3  Picobirnaviridae         Lysoka partiti-like virus     699     5
## 6      Iridoviridae      unclassified Chloriridovirus     575     6
## 4 Bacilladnaviridae     Marine protobacilladnavirus 1     503     7
## 1      unclassified  Guadeloupe Culex tymo-like virus 4168168     8
## 8      unclassified       Wenzhou sobemo-like virus 3 1066803     8
species_ord <- as.character(OTU_species_m_agg$species)
species_ord
## [1] "West Nile virus"                   "Phasi Charoen-like phasivirus"    
## [3] "Lampyris noctiluca errantivirus 1" "unclassified Circoviridae"        
## [5] "Lysoka partiti-like virus"         "unclassified Chloriridovirus"     
## [7] "Marine protobacilladnavirus 1"     "Guadeloupe Culex tymo-like virus" 
## [9] "Wenzhou sobemo-like virus 3"
which(species_ord == "West Nile virus")
## [1] 1
species_ord_1 <- c("West Nile virus", species_ord[-1])
species_ord_1
## [1] "West Nile virus"                   "Phasi Charoen-like phasivirus"    
## [3] "Lampyris noctiluca errantivirus 1" "unclassified Circoviridae"        
## [5] "Lysoka partiti-like virus"         "unclassified Chloriridovirus"     
## [7] "Marine protobacilladnavirus 1"     "Guadeloupe Culex tymo-like virus" 
## [9] "Wenzhou sobemo-like virus 3"
#### transform the table - reads for each viral species ##### 
rm <- c("Row.names", "superkingdom", "phylum", "class", "order", "family", "genus", "species")
OTU_species <- OTU_species[,!colnames(OTU_species) %in% rm]

#log
OTU_species_trans <- log10(OTU_species)
is.na(OTU_species_trans) <- sapply(OTU_species_trans, is.infinite)
OTU_species_trans[is.na(OTU_species_trans)]<-0
df.original <- OTU_species
df <- OTU_species_trans

#### use meta data for top annotation #####
meta = read.csv("~/Desktop/mosq_Guad_infection/GIM_metadata_new.csv", header=TRUE, row.names=1, sep=",")
meta_s = meta[!meta$mosq_species == "",]
meta_CQ <- meta_s[colnames(df),]

#### order sample by WNV reads within each group & column annotation  #### 
meta_CQz <- merge(meta_CQ, as.data.frame(t(df)),by=0)
meta_CQz$qPCR_body._ZIKV_WNV <- log10(meta_CQz$qPCR_body._ZIKV_WNV)

which(colnames(meta_CQz)=="qPCR_body._ZIKV_WNV")
## [1] 17
is.na(meta_CQz[,17])<-sapply(meta_CQz[,17], is.infinite)
meta_CQz[,17][is.na(meta_CQz[,17])]<-0
meta_CQz[,17]
##  [1] 4.209086 5.472123 3.452553 3.541704 4.729788 3.716671 4.211761 0.000000
##  [9] 2.972203 4.242566 3.451479 0.000000 3.004321 0.000000 2.136721 4.681802
## [17] 0.000000 3.691258 3.618048 4.242566 4.335879 0.000000 0.000000 4.521177
## [25] 4.574945 0.000000 4.350326 5.005627 0.000000 4.113609 3.787744 0.000000
## [33] 0.000000 4.256092 2.775246 4.221336 0.000000 0.000000 0.000000 0.000000
## [41] 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## [49] 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
meta_CQz$plaque_assay
##  [1] "7dpi-WNV(-)"      "7dpi-WNV(-)"      "7dpi-WNV(-)"      "7dpi-WNV(+)"     
##  [5] "7dpi-WNV(+)"      "7dpi-WNV(-)"      "7dpi-WNV(-)"      "7dpi-WNV(+)"     
##  [9] "7dpi-WNV(-)"      "7dpi-WNV(-)"      "7dpi-WNV(-)"      "7dpi-WNV(-)"     
## [13] "7dpi-WNV(-)"      "7dpi-WNV(-)"      "7dpi-WNV(-)"      "7dpi-WNV(-)"     
## [17] "14dpi-WNV(+)"     "14dpi-WNV(-)"     "14dpi-WNV(-)"     "14dpi-WNV(-)"    
## [21] "14dpi-WNV(-)"     "14dpi-WNV(-)"     "14dpi-WNV(-)"     "14dpi-WNV(+)"    
## [25] "14dpi-WNV(+)"     "14dpi-WNV(-)"     "14dpi-WNV(-)"     "14dpi-WNV(+)"    
## [29] "14dpi-WNV(-)"     "14dpi-WNV(-)"     "14dpi-WNV(+)"     "14dpi-WNV(-)"    
## [33] "14dpi-WNV(-)"     "14dpi-WNV(-)"     "14dpi-WNV(-)"     "14dpi-WNV(-)"    
## [37] "CQ-7dpi-blood"    "CQ-7dpi-blood"    "CQ-7dpi-blood"    "CQ-7dpi-blood"   
## [41] "CQ-7dpi-blood"    "CQ-14dpi-blood"   "CQ-14dpi-blood"   "CQ-14dpi-blood"  
## [45] "CQ-14dpi-blood"   "CQ-14dpi-blood"   "CQ-7dpi-sucrose"  "CQ-7dpi-sucrose" 
## [49] "CQ-7dpi-sucrose"  "CQ-7dpi-sucrose"  "CQ-14dpi-sucrose" "CQ-14dpi-sucrose"
## [53] "CQ-14dpi-sucrose" "CQ-14dpi-sucrose" "CQ-14dpi-sucrose"
meta_CQz$plaque_assay <- factor(meta_CQz$plaque_assay,
                                levels = c("CQ-7dpi-sucrose", "CQ-7dpi-blood", "7dpi-WNV(-)", "7dpi-WNV(+)", "CQ-14dpi-sucrose", "CQ-14dpi-blood", "14dpi-WNV(-)", "14dpi-WNV(+)"))

meta_CQz_order1 <- meta_CQz[meta_CQz$dpi== "7dpi",][with(meta_CQz[meta_CQz$dpi== "7dpi",], 
                                                         order(plaque_assay, -qPCR_body._ZIKV_WNV, -`West Nile virus`, -`Guadeloupe Culex tymo-like virus`,
                                                               -`Wenzhou sobemo-like virus 3`, -`Phasi Charoen-like phasivirus`)), ]
meta_CQz_order2 <- meta_CQz[meta_CQz$dpi== "14dpi",][with(meta_CQz[meta_CQz$dpi== "14dpi",], 
                                                          order(plaque_assay, -qPCR_body._ZIKV_WNV, -`West Nile virus`, -`Guadeloupe Culex tymo-like virus`,
                                                                -`Wenzhou sobemo-like virus 3`, -`Phasi Charoen-like phasivirus`)), ]
meta_CQz_order <- rbind(meta_CQz_order1,meta_CQz_order2)

meta_CQz$Row.names
##  [1] "GIM382" "GIM383" "GIM384" "GIM385" "GIM387" "GIM388" "GIM390" "GIM392"
##  [9] "GIM393" "GIM394" "GIM395" "GIM396" "GIM397" "GIM398" "GIM399" "GIM400"
## [17] "GIM421" "GIM422" "GIM423" "GIM424" "GIM425" "GIM426" "GIM427" "GIM428"
## [25] "GIM429" "GIM430" "GIM431" "GIM432" "GIM433" "GIM434" "GIM435" "GIM436"
## [33] "GIM437" "GIM438" "GIM439" "GIM440" "GIM451" "GIM452" "GIM453" "GIM454"
## [41] "GIM455" "GIM471" "GIM472" "GIM473" "GIM475" "GIM476" "GIM492" "GIM493"
## [49] "GIM494" "GIM495" "GIM525" "GIM527" "GIM528" "GIM529" "GIM530"
meta_CQz_order$Row.names
##  [1] "GIM495" "GIM493" "GIM494" "GIM492" "GIM452" "GIM451" "GIM453" "GIM454"
##  [9] "GIM455" "GIM383" "GIM400" "GIM394" "GIM390" "GIM382" "GIM388" "GIM384"
## [17] "GIM395" "GIM397" "GIM393" "GIM399" "GIM398" "GIM396" "GIM387" "GIM385"
## [25] "GIM392" "GIM525" "GIM527" "GIM529" "GIM530" "GIM528" "GIM471" "GIM475"
## [33] "GIM473" "GIM472" "GIM476" "GIM431" "GIM425" "GIM438" "GIM424" "GIM440"
## [41] "GIM434" "GIM422" "GIM423" "GIM439" "GIM426" "GIM436" "GIM433" "GIM430"
## [49] "GIM427" "GIM437" "GIM432" "GIM429" "GIM428" "GIM435" "GIM421"
df <- df[,meta_CQz_order$Row.names]
col_ha = HeatmapAnnotation(qPCR = anno_lines(meta_CQz_order[,17], gp = gpar(col = c("darkgrey")), 
                                             height = unit(1.5, "cm"),ylim = c(-0.5,6), axis_param = list(at = c(0,2,4,6)),
                                             add_points = TRUE, pt_gp = gpar(col = c("darkgrey")), pch = c(17)))

meta_CQz_order$group_n1 
## NULL
meta_CQz_order$plaque_assay
##  [1] CQ-7dpi-sucrose  CQ-7dpi-sucrose  CQ-7dpi-sucrose  CQ-7dpi-sucrose 
##  [5] CQ-7dpi-blood    CQ-7dpi-blood    CQ-7dpi-blood    CQ-7dpi-blood   
##  [9] CQ-7dpi-blood    7dpi-WNV(-)      7dpi-WNV(-)      7dpi-WNV(-)     
## [13] 7dpi-WNV(-)      7dpi-WNV(-)      7dpi-WNV(-)      7dpi-WNV(-)     
## [17] 7dpi-WNV(-)      7dpi-WNV(-)      7dpi-WNV(-)      7dpi-WNV(-)     
## [21] 7dpi-WNV(-)      7dpi-WNV(-)      7dpi-WNV(+)      7dpi-WNV(+)     
## [25] 7dpi-WNV(+)      CQ-14dpi-sucrose CQ-14dpi-sucrose CQ-14dpi-sucrose
## [29] CQ-14dpi-sucrose CQ-14dpi-sucrose CQ-14dpi-blood   CQ-14dpi-blood  
## [33] CQ-14dpi-blood   CQ-14dpi-blood   CQ-14dpi-blood   14dpi-WNV(-)    
## [37] 14dpi-WNV(-)     14dpi-WNV(-)     14dpi-WNV(-)     14dpi-WNV(-)    
## [41] 14dpi-WNV(-)     14dpi-WNV(-)     14dpi-WNV(-)     14dpi-WNV(-)    
## [45] 14dpi-WNV(-)     14dpi-WNV(-)     14dpi-WNV(-)     14dpi-WNV(-)    
## [49] 14dpi-WNV(-)     14dpi-WNV(-)     14dpi-WNV(+)     14dpi-WNV(+)    
## [53] 14dpi-WNV(+)     14dpi-WNV(+)     14dpi-WNV(+)    
## 8 Levels: CQ-7dpi-sucrose CQ-7dpi-blood 7dpi-WNV(-) ... 14dpi-WNV(+)
#### column split ####
df_1 <- df
df_1$WNV <- "others"
df_1$WNV[which(rownames(df_1)=="West Nile virus")] <- "aWNV"
#### draw figure #####
heatmap_col_tp2 <- colorRampPalette(brewer.pal(9, "YlGnBu"), space = "rgb")(100)
Heatmap(as.matrix(df),
        name = "log10 reads number", #title of legend
        column_title = "Samples", row_title = "viral species",
        row_names_gp = gpar(fontsize = 10),
        column_names_gp = gpar(fontsize = 6),
        column_title_gp = gpar(fontsize = 12),
        cluster_rows = F,
        cluster_columns= F,
        row_order = species_ord_1,
        column_labels = meta_CQz_order$mosq,
        top_annotation = col_ha,
        col = heatmap_col_tp2,
        column_split = meta_CQz_order$plaque_assay,
        row_split = df_1$WNV,
        border = TRUE,
        rect_gp = gpar(col = "white", lwd = 0.5)) 

2 qRT-PCR analysis

copies = read.table("~/Desktop/mosq_Guad_infection/qRT-PCR/qRTPCR_copies.csv", header=TRUE, dec=".", sep=",")
colnames(copies)
##  [1] "mosq_No"                           "mosq_species"                     
##  [3] "Treatment"                         "dpi"                              
##  [5] "part"                              "Body"                             
##  [7] "Head"                              "PCLPV_copies"                     
##  [9] "GMV_copies"                        "ATV_copies"                       
## [11] "AANV_copies"                       "ZIKV_copies"                      
## [13] "GCTLV_copies"                      "WSLV3_copies"                     
## [15] "WNV_copies"                        "group.plaque.qPCR.1000."          
## [17] "group.plaque.qPCR.0."              "group.plaque.qPCR.1000.body.1000."
## [19] "group.plaque.qPCR.0.body.0."       "Info1"                            
## [21] "plaque_assay"                      "head_ZIKV.1000"                   
## [23] "body_ZIKV.1000"
#### log transform ####
copies_log <- copies
copies_log <- cbind(copies[,1:7], log10(copies[,8:15]),copies[,16:23])
is.na(copies_log) <- sapply(copies_log, is.infinite)
copies_log[is.na(copies_log)] <- 0
colnames(copies_log)
##  [1] "mosq_No"                           "mosq_species"                     
##  [3] "Treatment"                         "dpi"                              
##  [5] "part"                              "Body"                             
##  [7] "Head"                              "PCLPV_copies"                     
##  [9] "GMV_copies"                        "ATV_copies"                       
## [11] "AANV_copies"                       "ZIKV_copies"                      
## [13] "GCTLV_copies"                      "WSLV3_copies"                     
## [15] "WNV_copies"                        "group.plaque.qPCR.1000."          
## [17] "group.plaque.qPCR.0."              "group.plaque.qPCR.1000.body.1000."
## [19] "group.plaque.qPCR.0.body.0."       "Info1"                            
## [21] "plaque_assay"                      "head_ZIKV.1000"                   
## [23] "body_ZIKV.1000"
########### all core viruses #########
copies_log$plaque_assay <- gsub("CQ-14dpi-","CQ-",copies_log$plaque_assay,fixed = T)
copies_log$plaque_assay <- gsub("CQ-7dpi-","CQ-",copies_log$plaque_assay,fixed = T)
copies_log$plaque_assay <- gsub("14dpi-","CQ-",copies_log$plaque_assay,fixed = T)
copies_log$plaque_assay <- gsub("7dpi-","CQ-",copies_log$plaque_assay,fixed = T)
copies_log$plaque_assay <- factor(copies_log$plaque_assay, levels = c("AA-sucrose", "AA-blood", "AA-ZIKV(-)", "AA-ZIKV(+)", 
                                                                      "CQ-sucrose", "CQ-blood", "CQ-WNV(-)", "CQ-WNV(+)"))
levels(copies_log$plaque_assay)
## [1] "AA-sucrose" "AA-blood"   "AA-ZIKV(-)" "AA-ZIKV(+)" "CQ-sucrose"
## [6] "CQ-blood"   "CQ-WNV(-)"  "CQ-WNV(+)"
copies_log$PCLPV_copies[copies_log$mosq_species == "Cx_quinquefasciatus"] <- NA
copies_log$GMV_copies[copies_log$mosq_species == "Cx_quinquefasciatus"] <- NA
copies_log$ATV_copies[copies_log$mosq_species == "Cx_quinquefasciatus"] <- NA
copies_log$AANV_copies[copies_log$mosq_species == "Cx_quinquefasciatus"] <- NA
copies_log$GCTLV_copies[copies_log$mosq_species == "Ae_aegypti"] <- NA
copies_log$WSLV3_copies[copies_log$mosq_species == "Ae_aegypti"] <- NA

copies_log_m <- melt(copies_log)
## Using mosq_No, mosq_species, Treatment, dpi, part, Body, Head, group.plaque.qPCR.1000., group.plaque.qPCR.0., group.plaque.qPCR.1000.body.1000., group.plaque.qPCR.0.body.0., Info1, plaque_assay, head_ZIKV.1000, body_ZIKV.1000 as id variables
copies_log_m_sub <- copies_log_m[copies_log_m$variable %in% c("PCLPV_copies", "GMV_copies", "ATV_copies","AANV_copies","GCTLV_copies", "WSLV3_copies"),]
copies_log_m_sub <- na.omit(copies_log_m_sub)

copies_log_m_sub$variable <- factor(copies_log_m_sub$variable, levels=c("GMV_copies", "PCLPV_copies","ATV_copies","AANV_copies","GCTLV_copies", "WSLV3_copies"))
levels(copies_log_m_sub$variable)
## [1] "GMV_copies"   "PCLPV_copies" "ATV_copies"   "AANV_copies"  "GCTLV_copies"
## [6] "WSLV3_copies"
copies_log_m_sub$virus[copies_log_m_sub$variable == "GMV_copies"] <- "Guadeloupe mosquito virus"
copies_log_m_sub$virus[copies_log_m_sub$variable == "PCLPV_copies"] <- "Phasi Charoen−like phasivirus"
copies_log_m_sub$virus[copies_log_m_sub$variable == "ATV_copies"] <- "Aedes aegypti totivirus"
copies_log_m_sub$virus[copies_log_m_sub$variable == "AANV_copies"] <- "Aedes anphevirus"
copies_log_m_sub$virus[copies_log_m_sub$variable == "GCTLV_copies"] <- "Guadeloupe Culex tymo−like virus"
copies_log_m_sub$virus[copies_log_m_sub$variable == "WSLV3_copies"] <- "Wenzhou sobemo−like virus 3"
copies_log_m_sub$virus <- factor(copies_log_m_sub$virus, levels=c("Guadeloupe mosquito virus", "Phasi Charoen−like phasivirus","Aedes aegypti totivirus",
                                                                  "Aedes anphevirus","Guadeloupe Culex tymo−like virus", "Wenzhou sobemo−like virus 3"))

copies_log_m_sub$dpi <- gsub("dpi", "dpe", copies_log_m_sub$dpi)
copies_log_m_sub$dpi<- factor(copies_log_m_sub$dpi, levels = c("7dpe", "21dpe", "14dpe"))

copies_log_m_sub_rmATV21 <- copies_log_m_sub[!(copies_log_m_sub$dpi=="21dpe"&copies_log_m_sub$variable=="ATV_copies"), ] 
copies_log_m_sub_rmGCTLV7 <- copies_log_m_sub_rmATV21[!(copies_log_m_sub_rmATV21$dpi=="7dpe"&copies_log_m_sub_rmATV21$variable=="GCTLV_copies"), ] 
copies_log_m_sub_rmWSLV37 <- copies_log_m_sub_rmGCTLV7[!(copies_log_m_sub_rmGCTLV7$dpi=="7dpe"&copies_log_m_sub_rmGCTLV7$variable=="WSLV3_copies"), ] 

col = plasma(8,alpha = 1, begin=0.9, end = 0, direction = 1)
ggplot(copies_log_m_sub_rmWSLV37, aes(x=plaque_assay, y=value, col=plaque_assay)) +
  geom_boxplot(na.rm = TRUE,outlier.shape = NA) + geom_jitter(shape=19, position=position_jitter(0.25), size = 2.2, alpha=0.7)+
  theme_bw()+
  facet_wrap(~ virus+part, scales="free", ncol = 6)+
  scale_color_manual(values = col) +
  theme(legend.position = "right",
        axis.ticks.x = element_blank(),
        axis.title.y = element_text(size=10), axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        strip.text.x = element_text(size = 12, colour = "black"),
        plot.title = element_text(hjust = 0.5, size = 12))+
  ggtitle("Viral genome copies determined by qRT-PCR")+
  ylab("log10 copies")+
  scale_y_continuous(limits = c(-0.01,9), breaks=seq(0, 9, by = 1))

copies_log_m_sub_rmATV7 <- copies_log_m_sub[!(copies_log_m_sub$dpi=="7dpe"&copies_log_m_sub$variable=="ATV_copies"), ] 
copies_log_m_sub_rmGCTLV21 <- copies_log_m_sub_rmATV7[!(copies_log_m_sub_rmATV7$dpi=="14dpe"&copies_log_m_sub_rmATV7$variable=="GCTLV_copies"), ] 
copies_log_m_sub_rmWSLV321 <- copies_log_m_sub_rmGCTLV21[!(copies_log_m_sub_rmGCTLV21$dpi=="14dpe"&copies_log_m_sub_rmGCTLV21$variable=="WSLV3_copies"), ] 

copies_log_m_sub_ATV <- copies_log_m_sub[copies_log_m_sub$variable=="ATV_copies", ]
copies_log_m_sub_ATV$dpi <- "7dpe+21dpe"
copies_log_m_sub_GCTLV <- copies_log_m_sub[copies_log_m_sub$variable=="GCTLV_copies", ]
copies_log_m_sub_GCTLV$dpi <- "7dpe+14dpe"
copies_log_m_sub_WSLV3 <- copies_log_m_sub[copies_log_m_sub$variable=="WSLV3_copies", ]
copies_log_m_sub_WSLV3$dpi <- "7dpe+14dpe"

copies_log_m_sub_new <- rbind(copies_log_m_sub_rmWSLV321,copies_log_m_sub_ATV)
copies_log_m_sub_new <- rbind(copies_log_m_sub_new,copies_log_m_sub_GCTLV)
copies_log_m_sub_new <- rbind(copies_log_m_sub_new,copies_log_m_sub_WSLV3)

ggplot(copies_log_m_sub_new, aes(x=plaque_assay, y=value, col=plaque_assay)) +
  geom_boxplot(na.rm = TRUE,outlier.shape = NA) + geom_jitter(shape=19, position=position_jitter(0.25), size = 2.2, alpha=0.7)+
  theme_bw()+
  facet_wrap(~ virus+dpi+part, scales="free", ncol = 4)+
  scale_color_manual(values = col) +
  theme(legend.position = "right",
        axis.ticks.x = element_blank(),
        axis.title.y = element_text(size=10), axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        strip.text.x = element_text(size = 12, colour = "black"),
        plot.title = element_text(hjust = 0.5, size = 12))+
  ggtitle("Viral genome copies determined by qRT-PCR")+
  ylab("log10 copies")+
  scale_y_continuous(limits = c(-0.01,9), breaks=seq(0, 9, by = 1))

########## pairwise.wilcox.test #######
#### GMV Guadeloupe mosquito virus
pairwise.wilcox.test(copies_log$GMV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="body"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="body"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  copies_log$GMV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "body"] 
## 
##            AA-sucrose AA-blood
## AA-blood   0.04491    -       
## AA-ZIKV(-) 1.2e-05    0.00059 
## 
## P value adjustment method: BH
pairwise.wilcox.test(copies_log$GMV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="head"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="head"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  copies_log$GMV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "head"] 
## 
##            AA-sucrose AA-blood
## AA-blood   0.0720     -       
## AA-ZIKV(-) 0.0011     0.0011  
## 
## P value adjustment method: BH
pairwise.wilcox.test(copies_log$GMV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="body"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="body"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  copies_log$GMV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "body"] 
## 
##            AA-sucrose AA-blood AA-ZIKV(-)
## AA-blood   0.0507     -        -         
## AA-ZIKV(-) 0.6847     0.0022   -         
## AA-ZIKV(+) 0.2667     0.0225   0.2667    
## 
## P value adjustment method: BH
pairwise.wilcox.test(copies_log$GMV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="head"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="head"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  copies_log$GMV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "head"] 
## 
##            AA-sucrose AA-blood AA-ZIKV(-)
## AA-blood   0.0507     -        -         
## AA-ZIKV(-) 0.0764     0.0031   -         
## AA-ZIKV(+) 0.3714     0.0225   0.9538    
## 
## P value adjustment method: BH
#### PCLPV Phasi Charoen−like phasivirus 
pairwise.wilcox.test(copies_log$PCLPV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="body"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="body"], 
                     p.adjust.method = "BH")
## 
##  Pairwise comparisons using Wilcoxon rank sum exact test 
## 
## data:  copies_log$PCLPV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "body"] 
## 
##            AA-sucrose AA-blood
## AA-blood   1.00       -       
## AA-ZIKV(-) 0.44       0.44    
## 
## P value adjustment method: BH
pairwise.wilcox.test(copies_log$PCLPV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="head"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="head"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  copies_log$PCLPV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "head"] 
## 
##            AA-sucrose AA-blood
## AA-blood   0.98       -       
## AA-ZIKV(-) 0.98       0.98    
## 
## P value adjustment method: BH
pairwise.wilcox.test(copies_log$PCLPV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="body"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="body"], 
                     p.adjust.method = "BH")
## 
##  Pairwise comparisons using Wilcoxon rank sum exact test 
## 
## data:  copies_log$PCLPV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "body"] 
## 
##            AA-sucrose AA-blood AA-ZIKV(-)
## AA-blood   0.84       -        -         
## AA-ZIKV(-) 0.35       0.32     -         
## AA-ZIKV(+) 0.35       0.32     0.46      
## 
## P value adjustment method: BH
pairwise.wilcox.test(copies_log$PCLPV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="head"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="head"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  copies_log$PCLPV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "head"] 
## 
##            AA-sucrose AA-blood AA-ZIKV(-)
## AA-blood   0.94       -        -         
## AA-ZIKV(-) 0.94       0.94     -         
## AA-ZIKV(+) 0.94       0.94     0.94      
## 
## P value adjustment method: BH
#### ATV Aedes aegypti totivirus 
pairwise.wilcox.test(copies_log$ATV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="body"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="body"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  copies_log$ATV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "body"] 
## 
##            AA-sucrose AA-blood
## AA-blood   0.11       -       
## AA-ZIKV(-) 0.03       0.52    
## 
## P value adjustment method: BH
pairwise.wilcox.test(copies_log$ATV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="head"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="head"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  copies_log$ATV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "head"] 
## 
##            AA-sucrose AA-blood
## AA-blood   0.100      -       
## AA-ZIKV(-) 0.061      1.000   
## 
## P value adjustment method: BH
pairwise.wilcox.test(copies_log$ATV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="body"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="body"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  copies_log$ATV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "body"] 
## 
##            AA-sucrose AA-blood AA-ZIKV(-)
## AA-blood   0.93       -        -         
## AA-ZIKV(-) 0.93       0.93     -         
## AA-ZIKV(+) 0.93       0.93     0.93      
## 
## P value adjustment method: BH
pairwise.wilcox.test(copies_log$ATV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="head"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="head"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  copies_log$ATV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "head"] 
## 
##            AA-sucrose AA-blood AA-ZIKV(-)
## AA-blood   0.359      -        -         
## AA-ZIKV(-) 0.366      0.039    -         
## AA-ZIKV(+) 0.824      0.216    0.366     
## 
## P value adjustment method: BH
#### AANV Aedes anphevirus
pairwise.wilcox.test(copies_log$AANV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="body"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="body"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  copies_log$AANV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "body"] 
## 
##            AA-sucrose AA-blood
## AA-blood   0.84       -       
## AA-ZIKV(-) 0.84       0.84    
## 
## P value adjustment method: BH
pairwise.wilcox.test(copies_log$AANV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="head"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="head"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  copies_log$AANV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "head"] 
## 
##            AA-sucrose AA-blood
## AA-blood   1          -       
## AA-ZIKV(-) 1          1       
## 
## P value adjustment method: BH
pairwise.wilcox.test(copies_log$AANV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="body"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="body"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  copies_log$AANV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "body"] 
## 
##            AA-sucrose AA-blood AA-ZIKV(-)
## AA-blood   1.00       -        -         
## AA-ZIKV(-) 0.58       0.68     -         
## AA-ZIKV(+) 0.58       0.58     0.58      
## 
## P value adjustment method: BH
pairwise.wilcox.test(copies_log$AANV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="head"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="head"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  copies_log$AANV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "head"] 
## 
##            AA-sucrose AA-blood AA-ZIKV(-)
## AA-blood   1.0        -        -         
## AA-ZIKV(-) 0.9        0.9      -         
## AA-ZIKV(+) 0.9        0.9      0.9       
## 
## P value adjustment method: BH
#### GCTLV Guadeloupe Culex tymo−like virus
pairwise.wilcox.test(copies_log$GCTLV_copies[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="7dpi" & copies_log$part=="body"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="7dpi" & copies_log$part=="body"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum exact test 
## 
## data:  copies_log$GCTLV_copies[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "7dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "7dpi" & copies_log$part == "body"] 
## 
##           CQ-sucrose CQ-blood CQ-WNV(-)
## CQ-blood  0.62       -        -        
## CQ-WNV(-) 0.62       1.00     -        
## CQ-WNV(+) 1.00       0.82     0.62     
## 
## P value adjustment method: BH
pairwise.wilcox.test(copies_log$GCTLV_copies[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="7dpi" & copies_log$part=="head"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="7dpi" & copies_log$part=="head"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  copies_log$GCTLV_copies[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "7dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "7dpi" & copies_log$part == "head"] 
## 
##           CQ-sucrose CQ-blood CQ-WNV(-)
## CQ-blood  0.553      -        -        
## CQ-WNV(-) 0.264      0.097    -        
## CQ-WNV(+) 0.264      0.048    0.553    
## 
## P value adjustment method: BH
pairwise.wilcox.test(copies_log$GCTLV_copies[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="14dpi" & copies_log$part=="body"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="14dpi" & copies_log$part=="body"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum exact test 
## 
## data:  copies_log$GCTLV_copies[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "14dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "14dpi" & copies_log$part == "body"] 
## 
##           CQ-sucrose CQ-blood CQ-WNV(-)
## CQ-blood  1.00       -        -        
## CQ-WNV(-) 0.59       0.59     -        
## CQ-WNV(+) 1.00       1.00     0.59     
## 
## P value adjustment method: BH
pairwise.wilcox.test(copies_log$GCTLV_copies[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="14dpi" & copies_log$part=="head"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="14dpi" & copies_log$part=="head"], 
                     p.adjust.method = "BH")
## 
##  Pairwise comparisons using Wilcoxon rank sum exact test 
## 
## data:  copies_log$GCTLV_copies[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "14dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "14dpi" & copies_log$part == "head"] 
## 
##           CQ-sucrose CQ-blood CQ-WNV(-)
## CQ-blood  0.421      -        -        
## CQ-WNV(-) 0.011      0.011    -        
## CQ-WNV(+) 0.016      0.048    0.202    
## 
## P value adjustment method: BH
#### WSLV3 Wenzhou sobemo−like virus 3
pairwise.wilcox.test(copies_log$WSLV3_copies[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="7dpi" & copies_log$part=="body"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="7dpi" & copies_log$part=="body"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  copies_log$WSLV3_copies[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "7dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "7dpi" & copies_log$part == "body"] 
## 
##           CQ-sucrose CQ-blood CQ-WNV(-)
## CQ-blood  0.18       -        -        
## CQ-WNV(-) 0.68       0.18     -        
## CQ-WNV(+) 0.51       0.15     0.18     
## 
## P value adjustment method: BH
pairwise.wilcox.test(copies_log$WSLV3_copies[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="7dpi" & copies_log$part=="head"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="7dpi" & copies_log$part=="head"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  copies_log$WSLV3_copies[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "7dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "7dpi" & copies_log$part == "head"] 
## 
##           CQ-sucrose CQ-blood CQ-WNV(-)
## CQ-blood  0.90       -        -        
## CQ-WNV(-) 0.90       0.90     -        
## CQ-WNV(+) 0.85       0.85     0.85     
## 
## P value adjustment method: BH
pairwise.wilcox.test(copies_log$WSLV3_copies[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="14dpi" & copies_log$part=="body"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="14dpi" & copies_log$part=="body"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  copies_log$WSLV3_copies[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "14dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "14dpi" & copies_log$part == "body"] 
## 
##           CQ-sucrose CQ-blood CQ-WNV(-)
## CQ-blood  0.242      -        -        
## CQ-WNV(-) 0.032      0.794    -        
## CQ-WNV(+) 0.797      0.149    0.032    
## 
## P value adjustment method: BH
pairwise.wilcox.test(copies_log$WSLV3_copies[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="14dpi" & copies_log$part=="head"], 
                     copies_log$plaque_assay[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="14dpi" & copies_log$part=="head"], 
                     p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  copies_log$WSLV3_copies[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "14dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "14dpi" & copies_log$part == "head"] 
## 
##           CQ-sucrose CQ-blood CQ-WNV(-)
## CQ-blood  0.447      -        -        
## CQ-WNV(-) 0.064      0.230    -        
## CQ-WNV(+) 0.797      0.447    0.064    
## 
## P value adjustment method: BH

3 16S rRNA analysis

#### load phyloseq project ####
GIMBac_decont <- readRDS("~/Desktop/mosq_Guad_infection/16s_rRNA/GIM_Bacteriome_decont_SLV132.rds")
GIMBac_decont
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 503 taxa and 158 samples ]
## sample_data() Sample Data:       [ 158 samples by 20 sample variables ]
## tax_table()   Taxonomy Table:    [ 503 taxa by 6 taxonomic ranks ]
col_6_Ae <-c("#8D634A", "#CCC591", "#84A0A5", "#972D15", "#C27D38", "#1D769C") 
col_6_Cq <-c("#7A7358",  "#85A789", "#C6CDF7","#29211F", "#376138","#7294D4")

3.1 Bacteriome of Aedes aegypti

GIMBac.Aa <- subset_samples(GIMBac_decont, mosq_species == "Ae_aegypti")
GIMBac.Aa <- subset_taxa(GIMBac.Aa, taxa_sums(GIMBac.Aa) > 0 )
GIMBac.Aa
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 381 taxa and 94 samples ]
## sample_data() Sample Data:       [ 94 samples by 20 sample variables ]
## tax_table()   Taxonomy Table:    [ 381 taxa by 6 taxonomic ranks ]

3.1.1 Genera proportion per sample

#### collape on Genus #### 
GIMBac.Aa.Genus <- GIMBac.Aa %>% #raw counts
  tax_glom(taxrank = "Genus")

GIMBac.Aa.Genus.per <- GIMBac.Aa.Genus  %>% 
  transform_sample_counts(function(x) {100*(x/sum(x))} ) %>%  
  psmelt()

GIMBac.Aa.Genus.per$group.plaque.qPCR.1000.body.1000.[which(is.na(GIMBac.Aa.Genus.per$group.plaque.qPCR.1000.body.1000.))] <- "NC"
levels(GIMBac.Aa.Genus.per$plaque_assay)
## NULL
#### 把一些低丰度的anno改成"others" #### 
AaBac_agg <- aggregate(Abundance ~ mosq + plaque_assay + Genus, FUN=sum, data=GIMBac.Aa.Genus.per)
AaBac_wide <- reshape2::dcast(AaBac_agg, mosq + plaque_assay ~ Genus, value.var = 'Abundance')
AaBac_wide <- AaBac_wide[order(-AaBac_wide$g_Asaia),]
sample_order <- AaBac_wide$mosq

AaBac_name_10 <- colnames(AaBac_wide[,3:ncol(AaBac_wide)][,colSums(AaBac_wide[,3:ncol(AaBac_wide)])<10])

nc <- which(colnames(GIMBac.Aa.Genus.per) == "Genus")
nc
## [1] 29
GIMBac.Aa.Genus.per[GIMBac.Aa.Genus.per$Genus %in% AaBac_name_10, nc] <- "Others"

GIMBac.Aa.Genus.per$Genus <- as.factor(GIMBac.Aa.Genus.per$Genus)
levels(GIMBac.Aa.Genus.per$Genus)
##  [1] "g_Acidovorax"            "g_Acinetobacter"        
##  [3] "g_Asaia"                 "g_Bacillus"             
##  [5] "g_Brevundimonas"         "g_Chryseobacterium"     
##  [7] "g_Cupriavidus"           "g_Klebsiella"           
##  [9] "g_Leucobacter"           "g_Mycobacterium"        
## [11] "g_Novosphingobium"       "g_Pseudomonas"          
## [13] "g_Salmonella"            "g_Sphingobacterium"     
## [15] "g_Sphingomonas"          "Others"                 
## [17] "uc_f_Enterobacteriaceae"
###--- repeat ---###
AaBac_agg <- aggregate(Abundance ~ mosq + plaque_assay + Genus, FUN=sum, data=GIMBac.Aa.Genus.per)
AaBac_wide <- reshape2::dcast(AaBac_agg, mosq + plaque_assay ~ Genus, value.var = 'Abundance')
as.data.frame(sort(colSums(AaBac_wide[,3:ncol(AaBac_wide)]) , decreasing=T))
##                         sort(colSums(AaBac_wide[, 3:ncol(AaBac_wide)]), decreasing = T)
## g_Asaia                                                                      6866.75066
## g_Chryseobacterium                                                            784.51212
## g_Pseudomonas                                                                 509.97242
## g_Novosphingobium                                                             371.14382
## g_Acidovorax                                                                  250.81813
## g_Klebsiella                                                                  206.27917
## g_Acinetobacter                                                               102.45281
## Others                                                                         79.11557
## g_Salmonella                                                                   45.36808
## uc_f_Enterobacteriaceae                                                        44.29459
## g_Cupriavidus                                                                  30.49819
## g_Mycobacterium                                                                29.13343
## g_Sphingomonas                                                                 19.51334
## g_Bacillus                                                                     17.45504
## g_Leucobacter                                                                  16.68823
## g_Sphingobacterium                                                             14.56847
## g_Brevundimonas                                                                11.43594
order <- rev(rownames(as.data.frame(sort(colSums(AaBac_wide[,3:ncol(AaBac_wide)]) , decreasing=T))))
which(order == "Others" | order == "uc_f_Enterobacteriaceae")
## [1]  8 10
order.1 <- c("Others", "uc_f_Enterobacteriaceae", order[-which(order == "Others" | order == "uc_f_Enterobacteriaceae")])
order.1
##  [1] "Others"                  "uc_f_Enterobacteriaceae"
##  [3] "g_Brevundimonas"         "g_Sphingobacterium"     
##  [5] "g_Leucobacter"           "g_Bacillus"             
##  [7] "g_Sphingomonas"          "g_Mycobacterium"        
##  [9] "g_Cupriavidus"           "g_Salmonella"           
## [11] "g_Acinetobacter"         "g_Klebsiella"           
## [13] "g_Acidovorax"            "g_Novosphingobium"      
## [15] "g_Pseudomonas"           "g_Chryseobacterium"     
## [17] "g_Asaia"
head(AaBac_agg)
##          mosq     plaque_assay        Genus Abundance
## 1 AaB-21dpi-1   AA-21dpi-blood g_Acidovorax         0
## 2 AaB-21dpi-2   AA-21dpi-blood g_Acidovorax         0
## 3 AaB-21dpi-3   AA-21dpi-blood g_Acidovorax         0
## 4 AaB-21dpi-4   AA-21dpi-blood g_Acidovorax         0
## 5 AaB-21dpi-5   AA-21dpi-blood g_Acidovorax         0
## 6 AaS-21dpi-1 AA-21dpi-sucrose g_Acidovorax         0
AaBac_agg$Genus <- factor(AaBac_agg$Genus, levels = order.1)
AaBac_agg$plaque_assay <- factor(AaBac_agg$plaque_assay, 
                                 levels = c("AA-7dpi-sucrose", "AA-7dpi-blood", "AA-7dpi-ZIKV(-)","AA-21dpi-sucrose", "AA-21dpi-blood", "AA-21dpi-ZIKV(-)", "AA-21dpi-ZIKV(+)"))
AaBac_agg$mosq <- factor(AaBac_agg$mosq, levels = sample_order)

col <- c("#CBD588", "#5F7FC7", "orange","#DA5724", "#508578", "#CD9BCD",
         "#AD6F3B", "#673770","#D14285", "#652926", "#C84248", 
         "#8569D5", "#5E738F","#D1A33D", "#8A7C64", "#599861","grey")

ggplot(AaBac_agg, aes(x = mosq, y = Abundance, fill = Genus)) + 
  facet_grid(~ plaque_assay,scales="free",space = "free") +  #### Lay out panels in a grid
  geom_bar(stat = "identity")+ 
  scale_fill_manual(values = rev(col)) +
  theme_bw()+
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(size = 8, angle = 90, vjust = 0.5,hjust = 1),
        axis.title.y = element_text(size = 12,vjust = 0.5),
        plot.margin = unit(c(1, 1, 1, 1), "cm"),
        strip.text.x = element_text(size =8, colour = "black"),
        plot.title = element_text(hjust = 0.5))

3.1.2 Alpha diversity on genus level

#### rarefy counts ####
sort(rowSums(otu_table(GIMBac.Aa)))
## GIM-320 GIM-318 GIM-319 GIM-317 GIM-549  GIM-30 GIM-257 GIM-259 GIM-580 GIM-586 
##     787    1445    1653    3635    4538    4978    9910   11428   16039   19805 
##  GIM-34 GIM-543  GIM-37 GIM-553  GIM-31 GIM-628 GIM-614 GIM-609 GIM-350  GIM-99 
##   20068   20530   21983   28947   33256   36063   42868   44769   45373   47031 
##  GIM-97 GIM-572  GIM-98  GIM-92 GIM-575 GIM-616  GIM-86  GIM-21  GIM-94 GIM-100 
##   47309   48020   48458   48841   50160   50288   50498   51337   51395   51467 
## GIM-627 GIM-618 GIM-617  GIM-36 GIM-260 GIM-611  GIM-26  GIM-95  GIM-24  GIM-32 
##   51962   52368   52446   52954   53403   53839   54581   55520   55666   55936 
##  GIM-83  GIM-82  GIM-87 GIM-574 GIM-290  GIM-84  GIM-29 GIM-582 GIM-585  GIM-81 
##   56574   56680   57223   57323   57427   57896   58005   58291   58441   59844 
##  GIM-38 GIM-347 GIM-546  GIM-39 GIM-621 GIM-288  GIM-35 GIM-610  GIM-89 GIM-625 
##   59877   60555   60803   61113   61342   61460   62137   62399   62488   62557 
## GIM-287  GIM-23 GIM-578 GIM-615  GIM-33 GIM-612 GIM-622 GIM-623 GIM-573 GIM-620 
##   63010   63549   63662   63813   64146   65068   65498   65642   67201   67355 
## GIM-284 GIM-583 GIM-576 GIM-577 GIM-258  GIM-22 GIM-608 GIM-581  GIM-96  GIM-28 
##   67642   67652   68252   68459   68582   68584   68685   68986   69076   70193 
## GIM-349  GIM-88 GIM-626  GIM-25 GIM-554 GIM-619  GIM-93  GIM-27  GIM-85  GIM-91 
##   70253   71987   72692   72716   73162   73723   74380   74467   74633   76212 
## GIM-613  GIM-90  GIM-40 GIM-624 
##   76585   82853   83026   90098
set.seed(100)
GIMBac.Aa_rare = rarefy_even_depth(GIMBac.Aa, sample.size = 11428, replace=FALSE, trimOTUs = F)
## You set `rngseed` to FALSE. Make sure you've set & recorded
##  the random seed of your session for reproducibility.
## See `?set.seed`
## ...
## 7 samples removedbecause they contained fewer reads than `sample.size`.
## Up to first five removed samples are:
## GIM-318GIM-30GIM-319GIM-320GIM-549
## ...
GIMBac.Aa_rm = prune_samples(sample_sums(GIMBac.Aa)<11428, GIMBac.Aa)
GIMBac.Aa_rm
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 381 taxa and 7 samples ]
## sample_data() Sample Data:       [ 7 samples by 20 sample variables ]
## tax_table()   Taxonomy Table:    [ 381 taxa by 6 taxonomic ranks ]
GIMBac.Aa_rareall <- merge_phyloseq(GIMBac.Aa_rare, GIMBac.Aa_rm)
sort(rowSums(otu_table(GIMBac.Aa_rareall)))
## GIM-320 GIM-318 GIM-319 GIM-317 GIM-549  GIM-30 GIM-257  GIM-21  GIM-29  GIM-36 
##     787    1445    1653    3635    4538    4978    9910   11428   11428   11428 
## GIM-575 GIM-585  GIM-87  GIM-94 GIM-609 GIM-617 GIM-625 GIM-543  GIM-22  GIM-37 
##   11428   11428   11428   11428   11428   11428   11428   11428   11428   11428 
## GIM-576 GIM-586  GIM-88  GIM-95 GIM-610 GIM-618 GIM-626 GIM-284  GIM-23  GIM-38 
##   11428   11428   11428   11428   11428   11428   11428   11428   11428   11428 
## GIM-577  GIM-81  GIM-89  GIM-96 GIM-611 GIM-619 GIM-627 GIM-287  GIM-24  GIM-31 
##   11428   11428   11428   11428   11428   11428   11428   11428   11428   11428 
##  GIM-39 GIM-578  GIM-82  GIM-97 GIM-612 GIM-620 GIM-628 GIM-288  GIM-25  GIM-32 
##   11428   11428   11428   11428   11428   11428   11428   11428   11428   11428 
##  GIM-40 GIM-580  GIM-83  GIM-90  GIM-98 GIM-613 GIM-621 GIM-290 GIM-347  GIM-26 
##   11428   11428   11428   11428   11428   11428   11428   11428   11428   11428 
##  GIM-33 GIM-572 GIM-581  GIM-84  GIM-91  GIM-99 GIM-614 GIM-622 GIM-258 GIM-349 
##   11428   11428   11428   11428   11428   11428   11428   11428   11428   11428 
##  GIM-27  GIM-34 GIM-573 GIM-582  GIM-85  GIM-92 GIM-100 GIM-615 GIM-623 GIM-259 
##   11428   11428   11428   11428   11428   11428   11428   11428   11428   11428 
## GIM-546 GIM-350  GIM-28  GIM-35 GIM-574 GIM-583  GIM-86  GIM-93 GIM-608 GIM-616 
##   11428   11428   11428   11428   11428   11428   11428   11428   11428   11428 
## GIM-624 GIM-260 GIM-553 GIM-554 
##   11428   11428   11428   11428
#### collape on Genus #### 
GIMBac.Aa_rare.Genus <- GIMBac.Aa_rareall %>% #raw counts
  tax_glom(taxrank = "Genus")

#### alpha diversity plot ####
plot_richness <- function (physeq, x = "samples", color = NULL, shape = NULL,
                           title = NULL, scales = "free_y", nrow = 1, shsi = NULL, measures = NULL,
                           sortby = NULL)
{
  erDF = estimate_richness(physeq, split = TRUE, measures = measures)
  measures = colnames(erDF)
  ses = colnames(erDF)[grep("^se\\.", colnames(erDF))]
  measures = measures[!measures %in% ses]
  if (!is.null(sample_data(physeq, errorIfNULL = FALSE))) {
    DF <- data.frame(erDF, sample_data(physeq))
  }
  else {
    DF <- data.frame(erDF)
  }
  if (!"samples" %in% colnames(DF)) {
    DF$samples <- sample_names(physeq)
  }
  if (!is.null(x)) {
    if (x %in% c("sample", "samples", "sample_names", "sample.names")) {
      x <- "samples"
    }
  }
  else {
    x <- "samples"
  }
  mdf = reshape2::melt(DF, measure.vars = measures)
  mdf$se <- NA_integer_
  if (length(ses) > 0) {
    selabs = ses
    names(selabs) <- substr(selabs, 4, 100)
    substr(names(selabs), 1, 1) <- toupper(substr(names(selabs),
                                                  1, 1))
    mdf$wse <- sapply(as.character(mdf$variable), function(i,
                                                           selabs) {
      selabs[i]
    }, selabs)
    for (i in 1:nrow(mdf)) {
      if (!is.na(mdf[i, "wse"])) {
        mdf[i, "se"] <- mdf[i, (mdf[i, "wse"])]
      }
    }
    mdf <- mdf[, -which(colnames(mdf) %in% c(selabs, "wse"))]
  }
  if (!is.null(measures)) {
    if (any(measures %in% as.character(mdf$variable))) {
      mdf <- mdf[as.character(mdf$variable) %in% measures,
      ]
    }
    else {
      warning("Argument to `measures` not supported. All alpha-diversity measures (should be) included in plot.")
    }
  }
  if (!is.null(shsi)) {
    warning("shsi no longer supported option in plot_richness. Please use `measures` instead")
  }
  if (!is.null(sortby)) {
    if (!all(sortby %in% levels(mdf$variable))) {
      warning("`sortby` argument not among `measures`. Ignored.")
    }
    if (!is.discrete(mdf[, x])) {
      warning("`sortby` argument provided, but `x` not a discrete variable. `sortby` is ignored.")
    }
    if (all(sortby %in% levels(mdf$variable)) & is.discrete(mdf[,
                                                                x])) {
      wh.sortby = which(mdf$variable %in% sortby)
      mdf[, x] <- factor(mdf[, x], levels = names(sort(tapply(X = mdf[wh.sortby,
                                                                      "value"], INDEX = mdf[wh.sortby, x], mean, na.rm = TRUE,
                                                              simplify = TRUE))))
    }
  }
  richness_map = aes_string(x = x, y = "value", colour = color,
                            shape = shape)
  p = ggplot(mdf, richness_map) + geom_boxplot(na.rm = TRUE,outlier.shape = NA)+geom_jitter(shape=19, position=position_jitter(0.25), size = 2.2, alpha=0.7)
  # if (any(!is.na(mdf[, "se"]))) {
  #   p = p + geom_errorbar(aes(ymax = value + se, ymin = value -
  #                               se), width = 0.1)
  # }
  p = p + theme(axis.text.x = element_text(angle = -90, vjust = 0.5,
                                           hjust = 0))
  p = p + ylab("Alpha Diversity Measure")
  p = p + facet_wrap(~variable, nrow = nrow, scales = scales)
  if (!is.null(title)) {
    p <- p + ggtitle(title)
  }
  return(p)
}

sample_data(GIMBac.Aa_rare.Genus)$Info1 <- factor(sample_data(GIMBac.Aa_rare.Genus)$Info1, levels = c("AA-7dpi-sucrose", "AA-7dpi-blood", "AA-7dpi-ZIKV", 
                                                                                                      "AA-21dpi-sucrose","AA-21dpi-blood", "AA-21dpi-ZIKV"))
sample_data(GIMBac.Aa_rare.Genus)$Treatment <- factor(sample_data(GIMBac.Aa_rare.Genus)$Treatment, levels = c("AA-sucrose", "AA-blood", "AA-ZIKV"))
sample_data(GIMBac.Aa_rare.Genus)$dpi <- factor(sample_data(GIMBac.Aa_rare.Genus)$dpi, levels = c("7dpi", "21dpi"))
GIMBac.Aa_rare.Genus.7d = subset_samples(GIMBac.Aa_rare.Genus, dpi == "7dpi")
GIMBac.Aa_rare.Genus.21d = subset_samples(GIMBac.Aa_rare.Genus, dpi == "21dpi")

plot_richness(GIMBac.Aa_rare.Genus, x="Info1",measures=c("Shannon"),col="Treatment", shape = "dpi") +
  scale_color_manual(values = rev(col_6_Ae[c(4:6)])) +
  theme_bw() +
  ylim(-0.001,1.9)+
  theme(axis.title.x = element_blank(),
        axis.text.x  = element_text(size = 12, colour = "black", angle = 90, hjust = 1, vjust=0.5),
        strip.text.x = element_text(size = 12, colour = "black"),
        plot.title = element_text(hjust = 0.5, size = 12)) +
  ggtitle("Alpha diversity of bacteriome in Ae. aegypti on genus level")+
  ylab("Alpha Diversity Value")

erich <- estimate_richness(GIMBac.Aa_rare.Genus, measures = c("Shannon"))
pairwise.wilcox.test(erich$Shannon,sample_data(GIMBac.Aa_rare.Genus)$Info1, p.adjust.method = "BH")
## 
##  Pairwise comparisons using Wilcoxon rank sum exact test 
## 
## data:  erich$Shannon and sample_data(GIMBac.Aa_rare.Genus)$Info1 
## 
##                  AA-7dpi-sucrose AA-7dpi-blood AA-7dpi-ZIKV AA-21dpi-sucrose
## AA-7dpi-blood    0.0433          -             -            -               
## AA-7dpi-ZIKV     0.2523          0.0075        -            -               
## AA-21dpi-sucrose 0.0433          0.5867        0.0433       -               
## AA-21dpi-blood   0.0170          0.0298        3.0e-05      0.0170          
## AA-21dpi-ZIKV    0.0028          0.9455        2.3e-08      0.2916          
##                  AA-21dpi-blood
## AA-7dpi-blood    -             
## AA-7dpi-ZIKV     -             
## AA-21dpi-sucrose -             
## AA-21dpi-blood   -             
## AA-21dpi-ZIKV    0.0028        
## 
## P value adjustment method: BH

3.1.3 Beta diversity on genus level

GIMBac.Aa_rare.GenusS <- subset_samples(GIMBac.Aa_rare.Genus, group.plaque.qPCR.1000.body.1000. != "")
sample_data(GIMBac.Aa_rare.GenusS)$group.plaque.qPCR.1000.body.1000. <- factor(sample_data(GIMBac.Aa_rare.GenusS)$group.plaque.qPCR.1000.body.1000.,
                                                                               levels =c("AA-7dpi-ZIKV(-/+/+)", "AA-7dpi-ZIKV(-/-/+)", "AA-7dpi-ZIKV(-/-/-)", 
                                                                                         "AA-7dpi-blood", "AA-7dpi-sucrose", "AA-21dpi-ZIKV(+/+/+)", 
                                                                                         "AA-21dpi-ZIKV(-/+/+)", "AA-21dpi-ZIKV(-/-/+)", "AA-21dpi-ZIKV(-/-/-)",
                                                                                         "AA-21dpi-blood", "AA-21dpi-sucrose"))
GIMBac.Aa_rare.Genus.7dS = subset_samples(GIMBac.Aa_rare.GenusS, dpi == "7dpi")
GIMBac.Aa_rare.Genus.21dS = subset_samples(GIMBac.Aa_rare.GenusS, dpi == "21dpi")
#### 7dpi ##### 
levels(sample_data(GIMBac.Aa_rare.Genus.7dS)$group.plaque.qPCR.1000.body.1000.)
## [1] "AA-7dpi-ZIKV(-/+/+)" "AA-7dpi-ZIKV(-/-/+)" "AA-7dpi-ZIKV(-/-/-)"
## [4] "AA-7dpi-blood"       "AA-7dpi-sucrose"
levels(sample_data(GIMBac.Aa_rare.Genus.7dS)$Treatment)
## [1] "AA-sucrose" "AA-blood"   "AA-ZIKV"
GP.ord.GIMBac.Aa_rare.Genus.7d <- ordinate(GIMBac.Aa_rare.Genus.7dS, "PCoA", "bray")

plot_ordination(GIMBac.Aa_rare.Genus.7dS, GP.ord.GIMBac.Aa_rare.Genus.7d, type = "samples", color = "Treatment", shape = "group.plaque.qPCR.1000.body.1000.")+
  theme_bw()+
  geom_point(size = 4)  + 
  scale_color_manual(values = col_6_Ae[c(4:6)])+
  scale_shape_manual(values=rev(shap[2:6])) +
  theme(legend.title=element_blank(),
        # legend.text= element_text(),
        legend.position = "right",
        plot.title = element_text(size = 12,hjust = 0.5))+
  ggtitle("PCoA of bacteriome in Ae. aegypti at 7 dpe") +
  stat_ellipse(type = "norm", linetype = 2,level=0.6,aes(group = Treatment))

set.seed(1)
# Calculate bray curtis distance matrix
GIMBac.Aa_rare.Genus.7dS_bray <- phyloseq::distance(GIMBac.Aa_rare.Genus.7dS, method = "bray")
# make a data frame from the sample_data
GIMBac.Aa_rare.Genus.7dS_sampledf <- data.frame(sample_data(GIMBac.Aa_rare.Genus.7dS))
# Adonis test
adonis(GIMBac.Aa_rare.Genus.7dS_bray ~ Treatment, data = GIMBac.Aa_rare.Genus.7dS_sampledf)
## 
## Call:
## adonis(formula = GIMBac.Aa_rare.Genus.7dS_bray ~ Treatment, data = GIMBac.Aa_rare.Genus.7dS_sampledf) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Treatment  2    2.8805 1.44024  10.932 0.35923  0.001 ***
## Residuals 39    5.1381 0.13175         0.64077           
## Total     41    8.0186                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#### 21dpi ##### 
levels(sample_data(GIMBac.Aa_rare.Genus.21dS)$group.plaque.qPCR.1000.body.1000.)
## [1] "AA-21dpi-ZIKV(+/+/+)" "AA-21dpi-ZIKV(-/+/+)" "AA-21dpi-ZIKV(-/-/+)"
## [4] "AA-21dpi-ZIKV(-/-/-)" "AA-21dpi-blood"       "AA-21dpi-sucrose"
levels(sample_data(GIMBac.Aa_rare.Genus.21dS)$Treatment)
## [1] "AA-sucrose" "AA-blood"   "AA-ZIKV"
GP.ord.GIMBac.Aa_rare.Genus.21d <- ordinate(GIMBac.Aa_rare.Genus.21dS, "PCoA", "bray")

plot_ordination(GIMBac.Aa_rare.Genus.21dS, GP.ord.GIMBac.Aa_rare.Genus.21d, type = "samples", color = "Treatment", shape = "group.plaque.qPCR.1000.body.1000.")+
  theme_bw()+
  geom_point(size = 4)  + 
  scale_color_manual(values = col_6_Ae[c(4:6)])+
  scale_shape_manual(values= c(shap[1:4],shap[6],shap[5])) +
  theme(legend.title=element_blank(),
        legend.position = "right",
        plot.title = element_text(size = 12,hjust = 0.5))+
  ggtitle("PCoA of bacteriome in Ae. aegypti at 21 dpe")

set.seed(1)
# Calculate bray curtis distance matrix
GIMBac.Aa_rare.Genus.21dS_bray <- phyloseq::distance(GIMBac.Aa_rare.Genus.21dS, method = "bray")
# make a data frame from the sample_data
GIMBac.Aa_rare.Genus.21dS_sampledf <- data.frame(sample_data(GIMBac.Aa_rare.Genus.21dS))
# Adonis test
adonis(GIMBac.Aa_rare.Genus.21dS_bray ~ Treatment, data = GIMBac.Aa_rare.Genus.21dS_sampledf)
## 
## Call:
## adonis(formula = GIMBac.Aa_rare.Genus.21dS_bray ~ Treatment,      data = GIMBac.Aa_rare.Genus.21dS_sampledf) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## Treatment  2   0.68676 0.34338  8.9531 0.27588  0.002 **
## Residuals 47   1.80259 0.03835         0.72412          
## Total     49   2.48935                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.2 Bacteriome of Culex quinquefasciatus

GIMBac.Cq <- subset_samples(GIMBac_decont, mosq_species == "Cx_quinquefasciatus")
GIMBac.Cq <- subset_taxa(GIMBac.Cq, taxa_sums(GIMBac.Cq) > 0 )
GIMBac.Cq
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 154 taxa and 60 samples ]
## sample_data() Sample Data:       [ 60 samples by 20 sample variables ]
## tax_table()   Taxonomy Table:    [ 154 taxa by 6 taxonomic ranks ]
sample_data(GIMBac.Cq)$group.plaque.qPCR.0.body.0. <- factor(sample_data(GIMBac.Cq)$group.plaque.qPCR.0.body.0., 
                                                             levels = c("CQ-7dpi-sucrose","CQ-7dpi-blood","7dpi-WNV(-/-)", "7dpi-WNV(-/+)", "7dpi-WNV(+/)", 
                                                                        "CQ-14dpi-sucrose","CQ-14dpi-blood","14dpi-WNV(-/-)","14dpi-WNV(-/+)","14dpi-WNV(+/)"))
sample_data(GIMBac.Cq)$Info1 <- factor(sample_data(GIMBac.Cq)$Info1, 
                                       levels = c("CQ-7dpi-sucrose", "CQ-7dpi-blood", "CQ-7dpi-WNV", "CQ-14dpi-sucrose", "CQ-14dpi-blood", "CQ-14dpi-WNV"))
sample_data(GIMBac.Cq)$Treatment <- factor(sample_data(GIMBac.Cq)$Treatment, levels = c("CQ-sucrose", "CQ-blood", "CQ-WNV"))

3.2.1 Genera proportion per sample

#### collape on Genus #### 
GIMBac.Cq.Genus <- GIMBac.Cq %>% #raw counts
  tax_glom(taxrank = "Genus")

GIMBac.Cq.Genus.per <- GIMBac.Cq.Genus  %>% 
  transform_sample_counts(function(x) {100*(x/sum(x))} ) %>%  
  psmelt()

colnames(GIMBac.Cq.Genus.per)
##  [1] "OTU"                               "Sample"                           
##  [3] "Abundance"                         "mosq"                             
##  [5] "body"                              "mosq_or"                          
##  [7] "head"                              "Index"                            
##  [9] "mosq_species"                      "Treatment"                        
## [11] "dpi"                               "group.plaque.qPCR.1000."          
## [13] "group.plaque.qPCR.0."              "group.plaque.qPCR.1000.body.1000."
## [15] "group.plaque.qPCR.0.body.0."       "Info1"                            
## [17] "plaque_assay"                      "qPCR_head._ZIKV_WNV"              
## [19] "qPCR_body._ZIKV_WNV"               "head_ZIKV.1000"                   
## [21] "body_ZIKV.1000"                    "Sample_or_Control"                
## [23] "conc"                              "Kingdom"                          
## [25] "Phylum"                            "Class"                            
## [27] "Order"                             "Family"                           
## [29] "Genus"
levels(GIMBac.Cq.Genus.per$plaque_assay)
## NULL
GIMBac.Cq.Genus.per$plaque_assay <- factor(GIMBac.Cq.Genus.per$plaque_assay, 
                                           levels = c("CQ-7dpi-sucrose", "CQ-7dpi-blood", "7dpi-WNV(-)", "7dpi-WNV(+)", 
                                                      "CQ-14dpi-sucrose", "CQ-14dpi-blood", "14dpi-WNV(-)", "14dpi-WNV(+)"))
levels(GIMBac.Cq.Genus.per$plaque_assay)
## [1] "CQ-7dpi-sucrose"  "CQ-7dpi-blood"    "7dpi-WNV(-)"      "7dpi-WNV(+)"     
## [5] "CQ-14dpi-sucrose" "CQ-14dpi-blood"   "14dpi-WNV(-)"     "14dpi-WNV(+)"
#### 把一些低丰度的anno改成“others”  ####
CqBac_agg <- aggregate(Abundance ~ mosq + plaque_assay + Genus, FUN=sum, data=GIMBac.Cq.Genus.per)
CqBac_wide <- dcast(CqBac_agg, mosq + plaque_assay ~ Genus, value.var = 'Abundance')
CqBac_wide_ordered <- CqBac_wide[with(CqBac_wide, order(plaque_assay,-g_Wolbachia,-g_Serratia)), ]
sample_order <- CqBac_wide_ordered$mosq

CqBac_name_10 <- colnames(CqBac_wide[,3:ncol(CqBac_wide)][,colSums(CqBac_wide[,3:ncol(CqBac_wide)])<10])

nc <- which(colnames(GIMBac.Cq.Genus.per) == "Genus")
nc
## [1] 29
GIMBac.Cq.Genus.per[GIMBac.Cq.Genus.per$Genus %in% CqBac_name_10, nc] <- "Others"

GIMBac.Cq.Genus.per$Genus <- as.factor(GIMBac.Cq.Genus.per$Genus)
levels(GIMBac.Cq.Genus.per$Genus)
## [1] "g_Asaia"            "g_Delftia"          "g_Elizabethkingia" 
## [4] "g_Gluconobacter"    "g_Pseudomonas"      "g_Serratia"        
## [7] "g_Stenotrophomonas" "g_Wolbachia"        "Others"
###--- repeat ---### 
CqBac_agg <- aggregate(Abundance ~ mosq + plaque_assay + Genus, FUN=sum, data=GIMBac.Cq.Genus.per)
CqBac_wide <- dcast(CqBac_agg, mosq + plaque_assay ~ Genus, value.var = 'Abundance')
as.data.frame(sort(colSums(CqBac_wide[,3:ncol(CqBac_wide)]) , decreasing=T))
##                    sort(colSums(CqBac_wide[, 3:ncol(CqBac_wide)]), decreasing = T)
## g_Wolbachia                                                             2789.01416
## g_Serratia                                                              1649.54264
## g_Asaia                                                                 1244.71360
## g_Elizabethkingia                                                        134.91983
## g_Pseudomonas                                                             71.33818
## g_Gluconobacter                                                           65.51299
## g_Stenotrophomonas                                                        16.96808
## Others                                                                    15.88106
## g_Delftia                                                                 12.10946
order <- rev(rownames(as.data.frame(sort(colSums(CqBac_wide[,3:ncol(CqBac_wide)]) , decreasing=T))))
which(order == "Others")
## [1] 2
order.1 <- c("Others", order[-which(order == "Others")])
order.1
## [1] "Others"             "g_Delftia"          "g_Stenotrophomonas"
## [4] "g_Gluconobacter"    "g_Pseudomonas"      "g_Elizabethkingia" 
## [7] "g_Asaia"            "g_Serratia"         "g_Wolbachia"
CqBac_agg$mosq <- factor(CqBac_agg$mosq, levels = sample_order)
CqBac_agg$Genus <- factor(CqBac_agg$Genus, levels = order.1)
levels(CqBac_agg$plaque_assay)
## [1] "CQ-7dpi-sucrose"  "CQ-7dpi-blood"    "7dpi-WNV(-)"      "7dpi-WNV(+)"     
## [5] "CQ-14dpi-sucrose" "CQ-14dpi-blood"   "14dpi-WNV(-)"     "14dpi-WNV(+)"
#### plot Genera proportion bar plot ####
col <- c("#5F7FC7", "#AD6F3B", "#CBD588", "#673770", "orange", "#CD9BCD", "#652926", "#C84248", "grey")
ggplot(CqBac_agg, aes(x = mosq, y = Abundance, fill = Genus)) + 
  facet_grid(~ plaque_assay,scales="free",space = "free") +  #### Lay out panels in a grid
  geom_bar(stat = "identity") + #### stat = "identity" : the height of the bar will represent the value in a column of the data frame
  scale_fill_manual(values = rev(col)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(size = 8, angle = 90, vjust = 0.5,hjust = 1),
        axis.title.y = element_text(size = 12,vjust = 0.5),
        plot.margin = unit(c(1, 1, 1, 1), "cm"),
        strip.text.x = element_text(size =8, colour = "black"),
        plot.title = element_text(hjust = 0.5)) + 
  guides(fill = guide_legend(reverse = TRUE, keywidth = 1, keyheight = 1))

3.2.2 Proportion of Asaia and Wolbachia

CqBac_agg$dpi <- CqBac_agg$plaque_assay
CqBac_agg$dpi <- gsub("CQ-", "", CqBac_agg$dpi)

CqBac_agg_v <- CqBac_agg %>% separate(dpi, c("dpi", "food"), sep ="-"); head(CqBac_agg_v)
## Warning: Expected 2 pieces. Additional pieces discarded in 270 rows [11, 12, 13,
## 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 41, 42, 43, 44, 45, ...].
##         mosq    plaque_assay   Genus  Abundance  dpi    food
## 1 CqS-7dpi-1 CQ-7dpi-sucrose g_Asaia 15.4421010 7dpi sucrose
## 2 CqS-7dpi-2 CQ-7dpi-sucrose g_Asaia  6.1206321 7dpi sucrose
## 3 CqS-7dpi-3 CQ-7dpi-sucrose g_Asaia  1.2262223 7dpi sucrose
## 4 CqS-7dpi-4 CQ-7dpi-sucrose g_Asaia  3.4048961 7dpi sucrose
## 5 CqS-7dpi-5 CQ-7dpi-sucrose g_Asaia  1.9477277 7dpi sucrose
## 6 CqB-7dpi-1   CQ-7dpi-blood g_Asaia  0.2025679 7dpi   blood
CqBac_agg_v$Genus <- gsub("g_", "", CqBac_agg_v$Genus); head(CqBac_agg_v)
##         mosq    plaque_assay Genus  Abundance  dpi    food
## 1 CqS-7dpi-1 CQ-7dpi-sucrose Asaia 15.4421010 7dpi sucrose
## 2 CqS-7dpi-2 CQ-7dpi-sucrose Asaia  6.1206321 7dpi sucrose
## 3 CqS-7dpi-3 CQ-7dpi-sucrose Asaia  1.2262223 7dpi sucrose
## 4 CqS-7dpi-4 CQ-7dpi-sucrose Asaia  3.4048961 7dpi sucrose
## 5 CqS-7dpi-5 CQ-7dpi-sucrose Asaia  1.9477277 7dpi sucrose
## 6 CqB-7dpi-1   CQ-7dpi-blood Asaia  0.2025679 7dpi   blood
CqBac_agg_v <- CqBac_agg_v[CqBac_agg_v$Genus %in% c("Asaia", "Wolbachia"), ]
CqBac_agg_v$dpi <- factor(CqBac_agg_v$dpi, levels= c("7dpi", "14dpi"))

wilcox.test(CqBac_agg_v[CqBac_agg_v$Genus=="Asaia"&CqBac_agg_v$dpi=="7dpi",]$Abundance,
            CqBac_agg_v[CqBac_agg_v$Genus=="Asaia"&CqBac_agg_v$dpi=="14dpi",]$Abundance, 
            p.adjust.method = "BH")
## 
##  Wilcoxon rank sum exact test
## 
## data:  CqBac_agg_v[CqBac_agg_v$Genus == "Asaia" & CqBac_agg_v$dpi == "7dpi", ]$Abundance and CqBac_agg_v[CqBac_agg_v$Genus == "Asaia" & CqBac_agg_v$dpi == "14dpi", ]$Abundance
## W = 182, p-value = 4.063e-05
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(CqBac_agg_v[CqBac_agg_v$Genus=="Wolbachia"&CqBac_agg_v$dpi=="7dpi",]$Abundance,
            CqBac_agg_v[CqBac_agg_v$Genus=="Wolbachia"&CqBac_agg_v$dpi=="14dpi",]$Abundance, 
            p.adjust.method = "BH")
## 
##  Wilcoxon rank sum exact test
## 
## data:  CqBac_agg_v[CqBac_agg_v$Genus == "Wolbachia" & CqBac_agg_v$dpi == "7dpi", ]$Abundance and CqBac_agg_v[CqBac_agg_v$Genus == "Wolbachia" & CqBac_agg_v$dpi == "14dpi", ]$Abundance
## W = 614, p-value = 0.01487
## alternative hypothesis: true location shift is not equal to 0
CqBac_agg_v$Food_sources[CqBac_agg_v$food == "WNV(+)"] <- "WNV(plaque+)"
CqBac_agg_v$Food_sources[CqBac_agg_v$food == "WNV("] <- "WNV(plaque-)"
CqBac_agg_v$Food_sources[CqBac_agg_v$food == "sucrose"] <- "sucrose"
CqBac_agg_v$Food_sources[CqBac_agg_v$food == "blood"] <- "blood"
CqBac_agg_v$Food_sources <- factor(CqBac_agg_v$Food_sources, levels=c("sucrose", "blood", "WNV(plaque-)", "WNV(plaque+)"))

ggplot(CqBac_agg_v, aes(x=Genus, y=Abundance,fill=dpi)) + 
  geom_boxplot(na.rm = TRUE,outlier.shape = NA)+
  facet_wrap(~Food_sources, nrow = 1)+
  theme_bw()+
  scale_fill_manual(values = col_6_Ae[c(3,1)])+
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(size = 12, hjust = 0.5),
        axis.text.y = element_text(size = 12),
        axis.title.y = element_text(size = 12,vjust = 0.5),
        # plot.margin = unit(c(1, 1, 1, 1), "cm"),
        strip.text.x = element_text(size =12, colour = "black"),
        plot.title = element_text(hjust = 0.5))+
  ylab("Proportion")+
  stat_compare_means()

3.2.2 Alpha diversity on genus level

#### rarefy  ####
sort(rowSums(otu_table(GIMBac.Cq)))
## GIM-434 GIM-495 GIM-530 GIM-527 GIM-440 GIM-400 GIM-452 GIM-382 GIM-396 GIM-392 
##    2367    6007   13219   24168   24676   29224   34447   38493   38678   39933 
## GIM-491 GIM-493 GIM-435 GIM-472 GIM-529 GIM-436 GIM-388 GIM-525 GIM-394 GIM-387 
##   42151   44527   46282   47873   49431   51469   51694   52552   53008   53342 
## GIM-473 GIM-528 GIM-391 GIM-433 GIM-475 GIM-454 GIM-381 GIM-430 GIM-437 GIM-492 
##   53406   53990   54060   56789   57476   57756   58002   58040   58072   58409 
## GIM-397 GIM-425 GIM-432 GIM-453 GIM-399 GIM-384 GIM-393 GIM-383 GIM-389 GIM-427 
##   58899   59979   60071   60274   60339   60551   60802   61024   61027   61520 
## GIM-424 GIM-471 GIM-428 GIM-431 GIM-421 GIM-451 GIM-385 GIM-390 GIM-422 GIM-438 
##   62791   62813   62965   63245   63659   64176   64967   65183   65903   67316 
## GIM-476 GIM-426 GIM-429 GIM-395 GIM-423 GIM-398 GIM-455 GIM-494 GIM-386 GIM-439 
##   68186   68512   68932   70712   71919   73396   73752   75039   76302   77901
set.seed(100)
GIMBac.Cq_rare = rarefy_even_depth(GIMBac.Cq, sample.size = 24168, replace=FALSE, trimOTUs = F)
## You set `rngseed` to FALSE. Make sure you've set & recorded
##  the random seed of your session for reproducibility.
## See `?set.seed`
## ...
## 3 samples removedbecause they contained fewer reads than `sample.size`.
## Up to first five removed samples are:
## GIM-495GIM-434GIM-530
## ...
GIMBac.Cq_rm = prune_samples(sample_sums(GIMBac.Cq)<24168, GIMBac.Cq)
GIMBac.Cq_rm
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 154 taxa and 3 samples ]
## sample_data() Sample Data:       [ 3 samples by 20 sample variables ]
## tax_table()   Taxonomy Table:    [ 154 taxa by 6 taxonomic ranks ]
GIMBac.Cq_rareall <- merge_phyloseq(GIMBac.Cq_rare, GIMBac.Cq_rm)

#### collape on Genus #### 
GIMBac.Cq_rare.Genus <- GIMBac.Cq_rareall %>% #raw counts
  tax_glom(taxrank = "Genus")

#### alpha diversity plot ####
GIMBac.Cq_rare.Genus
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 99 taxa and 60 samples ]
## sample_data() Sample Data:       [ 60 samples by 20 sample variables ]
## tax_table()   Taxonomy Table:    [ 99 taxa by 6 taxonomic ranks ]
GIMBac.Cq_rare.Genus.7d = subset_samples(GIMBac.Cq_rare.Genus, dpi == "7dpi")
GIMBac.Cq_rare.Genus.14d = subset_samples(GIMBac.Cq_rare.Genus, dpi == "14dpi")
sample_data(GIMBac.Cq_rare.Genus)$dpi <- factor(sample_data(GIMBac.Cq_rare.Genus)$dpi, levels = c("7dpi", "14dpi"))

plot_richness(GIMBac.Cq_rare.Genus, x="Info1",measures=c("Shannon"),col="Treatment", shape = "dpi") +
  scale_color_manual(values = rev(col_6_Ae[c(4:6)])) +
  theme_bw() +
  theme(axis.title.x = element_blank(),
        axis.text.x  = element_text(size = 12, colour = "black", angle = 90, hjust = 1, vjust=0.5),
        strip.text.x = element_text(size = 12, colour = "black"),
        plot.title = element_text(hjust = 0.5, size = 12)) +
  ggtitle("Alpha diversity of bacteriome in Cx. quinquefasciatus on genus level")+
  ylab("Alpha Diversity Value")

erich <- estimate_richness(GIMBac.Cq_rare.Genus, measures = c("Shannon"))
pairwise.wilcox.test(erich$Shannon,sample_data(GIMBac.Cq_rare.Genus)$Info1, p.adjust.method = "BH")
## 
##  Pairwise comparisons using Wilcoxon rank sum exact test 
## 
## data:  erich$Shannon and sample_data(GIMBac.Cq_rare.Genus)$Info1 
## 
##                  CQ-7dpi-sucrose CQ-7dpi-blood CQ-7dpi-WNV CQ-14dpi-sucrose
## CQ-7dpi-blood    0.32            -             -           -               
## CQ-7dpi-WNV      0.32            0.89          -           -               
## CQ-14dpi-sucrose 0.86            0.70          0.32        -               
## CQ-14dpi-blood   0.58            0.90          0.78        0.32            
## CQ-14dpi-WNV     0.97            0.32          0.12        0.78            
##                  CQ-14dpi-blood
## CQ-7dpi-blood    -             
## CQ-7dpi-WNV      -             
## CQ-14dpi-sucrose -             
## CQ-14dpi-blood   -             
## CQ-14dpi-WNV     0.14          
## 
## P value adjustment method: BH

3.2.3 Beta diversity on genus level

GIMBac.Cq_rare.Genus
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 99 taxa and 60 samples ]
## sample_data() Sample Data:       [ 60 samples by 20 sample variables ]
## tax_table()   Taxonomy Table:    [ 99 taxa by 6 taxonomic ranks ]
GIMBac.Cq_rare.Genus.7d
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 99 taxa and 30 samples ]
## sample_data() Sample Data:       [ 30 samples by 20 sample variables ]
## tax_table()   Taxonomy Table:    [ 99 taxa by 6 taxonomic ranks ]
GIMBac.Cq_rare.Genus.14d
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 99 taxa and 30 samples ]
## sample_data() Sample Data:       [ 30 samples by 20 sample variables ]
## tax_table()   Taxonomy Table:    [ 99 taxa by 6 taxonomic ranks ]
#### 7dpi ##### 
GP.ord.GIMBac.Cq_rare.Genus.7d <- ordinate(GIMBac.Cq_rare.Genus.7d, "PCoA", "bray")

plot_ordination(GIMBac.Cq_rare.Genus.7d, GP.ord.GIMBac.Cq_rare.Genus.7d, type = "samples", color = "Treatment", shape = "group.plaque.qPCR.0.body.0.")+
  theme_bw()+
  geom_point(size = 4)  + 
  scale_color_manual(values = col_6_Ae[c(4:6)])+
  scale_shape_manual(values=c(shap[1:4],shap[6])) +
  theme(legend.title=element_blank(),
        legend.position = "right",
        plot.title = element_text(size = 12,hjust = 0.5))+
  ggtitle("PCoA of bacteriome in Cx. quinquefasciatus at 7dpi")

dev.off()
## null device 
##           1
set.seed(1)
# Calculate bray curtis distance matrix
GIMBac.Cq_rare.Genus.7d_bray <- phyloseq::distance(GIMBac.Cq_rare.Genus.7d, method = "bray")
# make a data frame from the sample_data
GIMBac.Cq_rare.Genus.7d_sampledf <- data.frame(sample_data(GIMBac.Cq_rare.Genus.7d))
# Adonis test
adonis(GIMBac.Cq_rare.Genus.7d_bray ~ Treatment, data = GIMBac.Cq_rare.Genus.7d_sampledf)
## 
## Call:
## adonis(formula = GIMBac.Cq_rare.Genus.7d_bray ~ Treatment, data = GIMBac.Cq_rare.Genus.7d_sampledf) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model     R2 Pr(>F)
## Treatment  2    0.5371 0.26853  1.4751 0.0985  0.184
## Residuals 27    4.9153 0.18205         0.9015       
## Total     29    5.4524                 1.0000
#### 14dpi ##### 
GP.ord.GIMBac.Cq_rare.Genus.14d <- ordinate(GIMBac.Cq_rare.Genus.14d, "PCoA", "bray")

plot_ordination(GIMBac.Cq_rare.Genus.14d, GP.ord.GIMBac.Cq_rare.Genus.14d, type = "samples", color = "Treatment", shape = "group.plaque.qPCR.0.body.0.")+
  theme_bw()+
  geom_point(size = 4)  + 
  scale_color_manual(values = col_6_Ae[c(4:6)])+
  scale_shape_manual(values=c(shap[1:4],shap[6])) +
  theme(legend.title=element_blank(),
        legend.position = "right",
        plot.title = element_text(size = 12,hjust = 0.5))+
  ggtitle("PCoA of bacteriome in Cx. quinquefasciatus at 14dpi")

set.seed(1)
# Calculate bray curtis distance matrix
GIMBac.Cq_rare.Genus.14d_bray <- phyloseq::distance(GIMBac.Cq_rare.Genus.14d, method = "bray")
# make a data frame from the sample_data
GIMBac.Cq_rare.Genus.14d_sampledf <- data.frame(sample_data(GIMBac.Cq_rare.Genus.14d))
# Adonis test
adonis(GIMBac.Cq_rare.Genus.14d_bray ~ Treatment, data = GIMBac.Cq_rare.Genus.14d_sampledf)
## 
## Call:
## adonis(formula = GIMBac.Cq_rare.Genus.14d_bray ~ Treatment, data = GIMBac.Cq_rare.Genus.14d_sampledf) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## Treatment  2    0.4523 0.22613  1.1465 0.07828  0.353
## Residuals 27    5.3254 0.19724         0.92172       
## Total     29    5.7777                 1.00000

4 DESeq2 determined differential bacteria

#### Aedes ####
df_otu <- as.data.frame(t(GIMBac.Aa_rare.Genus@otu_table))
df_tax <- as.data.frame(GIMBac.Aa_rare.Genus@tax_table)
df_meta <- data.frame(sample_data(GIMBac.Aa_rare.Genus))
df <- merge(df_otu, df_tax, by = "row.names", all=T)
#### Aedes-Diet & Aedes-Infectious ####
AAdiet <- readRDS("~/Desktop/mosq_Guad_infection/RDS/bacteria_DESeq2_AA_treatment.RDS")
AAdiet
## [[1]]
##          baseMean log2FoldChange    lfcSE       stat       pvalue         padj
## SVs5   591.326427     -28.084810 1.525607 -18.408936 1.113902e-75 6.460630e-74
## SVs6   298.034542     -27.219452 1.647940 -16.517262 2.756314e-61 7.993311e-60
## SVs4  1492.902368     -13.216161 1.265541 -10.443093 1.575906e-25 3.046752e-24
## SVs15   35.077201      -8.350071 1.310547  -6.371438 1.872640e-10 2.715327e-09
## SVs12   36.391970      -4.582519 1.002896  -4.569287 4.893853e-06 5.676869e-05
## SVs45    3.197046      -4.904141 1.229969  -3.987206 6.685584e-05 6.462731e-04
##                    Genus           cat
## SVs5   g_Novosphingobium AAbloodvsZIKV
## SVs6        g_Acidovorax AAbloodvsZIKV
## SVs4  g_Chryseobacterium AAbloodvsZIKV
## SVs15      g_Cupriavidus AAbloodvsZIKV
## SVs12    g_Acinetobacter AAbloodvsZIKV
## SVs45      g_Roseococcus AAbloodvsZIKV
## 
## [[2]]
##           baseMean log2FoldChange     lfcSE       stat       pvalue
## SVs5    495.298867     -23.741984 1.5667411 -15.153738 7.157719e-52
## SVs6    249.635674     -25.546161 1.6928788 -15.090366 1.874008e-51
## SVs9   2277.909223      17.218830 2.1226043   8.112124 4.974231e-16
## SVs4   1254.221186      -5.649115 1.1792208  -4.790549 1.663254e-06
## SVs15    29.380892      -6.098126 1.3126416  -4.645690 3.389420e-06
## SVs78     4.132256       6.416490 1.5232939   4.212247 2.528427e-05
## SVs8    344.495062       2.854758 0.7157120   3.988697 6.643732e-05
## SVs35     9.598872       2.891754 0.8544597   3.384307 7.135821e-04
## SVs102    6.156045      11.388979 3.7627032   3.026808 2.471511e-03
## SVs71     6.554007       5.420686 1.9483073   2.782254 5.398278e-03
##                padj              Genus             cat
## SVs5   2.863088e-50  g_Novosphingobium AAsucrosevsZIKV
## SVs6   3.748015e-50       g_Acidovorax AAsucrosevsZIKV
## SVs9   6.632307e-15       g_Klebsiella AAsucrosevsZIKV
## SVs4   1.663254e-05 g_Chryseobacterium AAsucrosevsZIKV
## SVs15  2.711536e-05      g_Cupriavidus AAsucrosevsZIKV
## SVs78  1.685618e-04     g_Sphingomonas AAsucrosevsZIKV
## SVs8   3.796418e-04      g_Pseudomonas AAsucrosevsZIKV
## SVs35  3.567911e-03         g_Bacillus AAsucrosevsZIKV
## SVs102 1.098449e-02 g_Sphingobacterium AAsucrosevsZIKV
## SVs71  2.159311e-02    g_Brevundimonas AAsucrosevsZIKV
## 
## [[3]]
##         baseMean log2FoldChange     lfcSE     stat       pvalue         padj
## SVs102  22.03441      24.131719 2.8839593 8.367566 5.882061e-17 4.176263e-15
## SVs9   228.68067       8.504409 1.5953683 5.330687 9.784177e-08 3.473383e-06
## SVs8   647.40055       4.005112 0.9860531 4.061761 4.870402e-05 1.152662e-03
## SVs71   19.20429       8.782719 2.2848029 3.843972 1.210586e-04 2.148791e-03
## SVs12   14.66698       5.430173 1.5982411 3.397593 6.798141e-04 9.653360e-03
##                     Genus              cat
## SVs102 g_Sphingobacterium AAsucrosevsblood
## SVs9         g_Klebsiella AAsucrosevsblood
## SVs8        g_Pseudomonas AAsucrosevsblood
## SVs71     g_Brevundimonas AAsucrosevsblood
## SVs12     g_Acinetobacter AAsucrosevsblood
AAinfec <- readRDS("~/Desktop/mosq_Guad_infection/RDS/Bacteria_DESeq2_AA_infectious.RDS")
AAinfec
## [[1]]
##       baseMean log2FoldChange    lfcSE     stat       pvalue         padj
## SVs12 122.6242        7.87385 1.298181 6.065296 1.317112e-09 1.040519e-07
##                 Genus              cat
## SVs12 g_Acinetobacter 21plaqueposvsneg
gp <- as.data.frame(sort(table(c(AAdiet[[1]]$Genus, AAdiet[[2]]$Genus, AAdiet[[3]]$Genus, AAinfec[[1]]$Genus)), decreasing = T))
gp <- gp[-c(8,10,12),]
order <- c("Chryseobacterium", "Novosphingobium", "Acidovorax", "Cupriavidus", "Roseococcus",
           "Acinetobacter", "Klebsiella", "Brevundimonas", "Sphingobacterium")

df_sp <- df[df$Genus %in% as.character(gp$Var1),]
df_sp <- df_sp[, colnames(df_sp) %in% c(colnames(df_otu), "Genus")]
df_spm <- melt(df_sp); head(df_spm)
## Using Genus as id variables
##                Genus variable value
## 1 g_Sphingobacterium   GIM-21     0
## 2    g_Acinetobacter   GIM-21    17
## 3      g_Cupriavidus   GIM-21    42
## 4 g_Chryseobacterium   GIM-21  2012
## 5      g_Roseococcus   GIM-21     0
## 6  g_Novosphingobium   GIM-21    79
df_meta$variable <- rownames(df_meta)
df_merge <- merge(df_spm, df_meta, by = "variable", all =T); head(df_merge)
##   variable              Genus value         mosq    body      mosq_or   head
## 1  GIM-100 g_Sphingobacterium     0 AaZ-21dpi-20 GIM-100 AMT-21dpi-20 GIM-80
## 2  GIM-100    g_Acinetobacter     2 AaZ-21dpi-20 GIM-100 AMT-21dpi-20 GIM-80
## 3  GIM-100      g_Cupriavidus     0 AaZ-21dpi-20 GIM-100 AMT-21dpi-20 GIM-80
## 4  GIM-100 g_Chryseobacterium  5965 AaZ-21dpi-20 GIM-100 AMT-21dpi-20 GIM-80
## 5  GIM-100      g_Roseococcus     0 AaZ-21dpi-20 GIM-100 AMT-21dpi-20 GIM-80
## 6  GIM-100  g_Novosphingobium     0 AaZ-21dpi-20 GIM-100 AMT-21dpi-20 GIM-80
##           Index mosq_species Treatment   dpi group.plaque.qPCR.1000.
## 1 AMT-21dpi-B20   Ae_aegypti   AA-ZIKV 21dpi      AA-21dpi-ZIKV(-/-)
## 2 AMT-21dpi-B20   Ae_aegypti   AA-ZIKV 21dpi      AA-21dpi-ZIKV(-/-)
## 3 AMT-21dpi-B20   Ae_aegypti   AA-ZIKV 21dpi      AA-21dpi-ZIKV(-/-)
## 4 AMT-21dpi-B20   Ae_aegypti   AA-ZIKV 21dpi      AA-21dpi-ZIKV(-/-)
## 5 AMT-21dpi-B20   Ae_aegypti   AA-ZIKV 21dpi      AA-21dpi-ZIKV(-/-)
## 6 AMT-21dpi-B20   Ae_aegypti   AA-ZIKV 21dpi      AA-21dpi-ZIKV(-/-)
##   group.plaque.qPCR.0. group.plaque.qPCR.1000.body.1000.
## 1   AA-21dpi-ZIKV(-/-)              AA-21dpi-ZIKV(-/-/-)
## 2   AA-21dpi-ZIKV(-/-)              AA-21dpi-ZIKV(-/-/-)
## 3   AA-21dpi-ZIKV(-/-)              AA-21dpi-ZIKV(-/-/-)
## 4   AA-21dpi-ZIKV(-/-)              AA-21dpi-ZIKV(-/-/-)
## 5   AA-21dpi-ZIKV(-/-)              AA-21dpi-ZIKV(-/-/-)
## 6   AA-21dpi-ZIKV(-/-)              AA-21dpi-ZIKV(-/-/-)
##   group.plaque.qPCR.0.body.0.         Info1     plaque_assay
## 1        AA-21dpi-ZIKV(-/-/-) AA-21dpi-ZIKV AA-21dpi-ZIKV(-)
## 2        AA-21dpi-ZIKV(-/-/-) AA-21dpi-ZIKV AA-21dpi-ZIKV(-)
## 3        AA-21dpi-ZIKV(-/-/-) AA-21dpi-ZIKV AA-21dpi-ZIKV(-)
## 4        AA-21dpi-ZIKV(-/-/-) AA-21dpi-ZIKV AA-21dpi-ZIKV(-)
## 5        AA-21dpi-ZIKV(-/-/-) AA-21dpi-ZIKV AA-21dpi-ZIKV(-)
## 6        AA-21dpi-ZIKV(-/-/-) AA-21dpi-ZIKV AA-21dpi-ZIKV(-)
##   qPCR_head._ZIKV_WNV qPCR_body._ZIKV_WNV   head_ZIKV.1000   body_ZIKV.1000
## 1                   0                   0 AA-21dpi-ZIKV(-) AA-21dpi-ZIKV(-)
## 2                   0                   0 AA-21dpi-ZIKV(-) AA-21dpi-ZIKV(-)
## 3                   0                   0 AA-21dpi-ZIKV(-) AA-21dpi-ZIKV(-)
## 4                   0                   0 AA-21dpi-ZIKV(-) AA-21dpi-ZIKV(-)
## 5                   0                   0 AA-21dpi-ZIKV(-) AA-21dpi-ZIKV(-)
## 6                   0                   0 AA-21dpi-ZIKV(-) AA-21dpi-ZIKV(-)
##   Sample_or_Control    conc
## 1              TRUE 11.0088
## 2              TRUE 11.0088
## 3              TRUE 11.0088
## 4              TRUE 11.0088
## 5              TRUE 11.0088
## 6              TRUE 11.0088
df_merge$log10reads <- log10(df_merge$value)
is.na(df_merge$log10reads) <- sapply(df_merge$log10reads, is.infinite)
df_merge$log10reads[is.na(df_merge$log10reads)] <- 0
df_merge$Genus <- gsub("g_", "", df_merge$Genus)
df_merge$Genus <- factor(df_merge$Genus, levels = order); levels(df_merge$Genus)
## [1] "Chryseobacterium" "Novosphingobium"  "Acidovorax"       "Cupriavidus"     
## [5] "Roseococcus"      "Acinetobacter"    "Klebsiella"       "Brevundimonas"   
## [9] "Sphingobacterium"
df_merge_aa <- df_merge
df_merge_aa <- df_merge_aa[,colnames(df_merge_aa) %in% c("Treatment", "Genus", "variable", "log10reads")]

#### Culex ####
df_otu <- as.data.frame(t(GIMBac.Cq_rare.Genus@otu_table))
df_tax <- as.data.frame(GIMBac.Cq_rare.Genus@tax_table)
df_meta <- data.frame(sample_data(GIMBac.Cq_rare.Genus))
df <- merge(df_otu, df_tax, by = "row.names", all=T)
#### Culex-Diet ####
CQdiet <- readRDS("~/Desktop/mosq_Guad_infection/RDS/Bacteria_DESeq2_CQ_treatment.RDS")
CQdiet
## [[1]]
##      baseMean log2FoldChange    lfcSE     stat       pvalue        padj
## SVs3  70296.5       4.432407 1.163301 3.810196 0.0001388567 0.008886832
##           Genus          cat
## SVs3 g_Serratia CQbloodvsWNV
## 
## [[2]]
##         baseMean log2FoldChange     lfcSE     stat       pvalue         padj
## SVs7  325.366291      30.000000 3.6211361 8.284693 1.184144e-16 9.591566e-15
## SVs85   1.855982      23.153094 3.6296383 6.378898 1.783664e-10 7.223838e-09
## SVs8  305.228196       3.813213 0.7360724 5.180487 2.213078e-07 5.975312e-06
## SVs92   9.928171       8.105716 2.5808287 3.140741 1.685207e-03 3.412545e-02
##                   Genus            cat
## SVs7  g_Elizabethkingia CQsucrosevsWNV
## SVs85      g_Legionella CQsucrosevsWNV
## SVs8      g_Pseudomonas CQsucrosevsWNV
## SVs92  g_Pseudonocardia CQsucrosevsWNV
## 
## [[3]]
##        baseMean log2FoldChange    lfcSE     stat       pvalue         padj
## SVs85  9.093585      23.018909 2.960427 7.775538 7.512755e-15 4.582781e-13
## SVs35 36.037711       6.273055 1.320633 4.750038 2.033788e-06 6.203053e-05
## SVs55 15.306501       7.452474 1.913401 3.894885 9.824549e-05 1.997658e-03
## SVs92 14.322496       7.325056 2.475187 2.959395 3.082438e-03 4.700719e-02
##                  Genus              cat
## SVs85     g_Legionella CQsucrosevsblood
## SVs35       g_Bacillus CQsucrosevsblood
## SVs55      g_Pelomonas CQsucrosevsblood
## SVs92 g_Pseudonocardia CQsucrosevsblood
gp <- as.data.frame(sort(table(c(CQdiet[[1]]$Genus, CQdiet[[2]]$Genus, CQdiet[[3]]$Genus)), decreasing = T))
gp
##                Var1 Freq
## 1      g_Legionella    2
## 2  g_Pseudonocardia    2
## 3        g_Bacillus    1
## 4 g_Elizabethkingia    1
## 5       g_Pelomonas    1
## 6     g_Pseudomonas    1
## 7        g_Serratia    1
order <- c("Serratia", "Pseudomonas", "Bacillus")

df_sp <- df[df$Genus %in% c("g_Serratia", "g_Pseudomonas", "g_Bacillus"),]
df_sp <- df_sp[, colnames(df_sp) %in% c(colnames(df_otu), "Genus")]
df_spm <- melt(df_sp); head(df_spm)
## Using Genus as id variables
##           Genus variable value
## 1    g_Serratia  GIM-387     0
## 2    g_Bacillus  GIM-387    52
## 3 g_Pseudomonas  GIM-387    99
## 4    g_Serratia  GIM-394  7252
## 5    g_Bacillus  GIM-394     6
## 6 g_Pseudomonas  GIM-394   100
df_meta$variable <- rownames(df_meta)
df_merge <- merge(df_spm, df_meta, by = "variable", all =T); head(df_merge)
##   variable         Genus value       mosq    body    mosq_or    head
## 1  GIM-381    g_Serratia 20620 CqW-7dpi-1 GIM-381 CWV-7dpi-1 GIM-361
## 2  GIM-381    g_Bacillus     5 CqW-7dpi-1 GIM-381 CWV-7dpi-1 GIM-361
## 3  GIM-381 g_Pseudomonas    14 CqW-7dpi-1 GIM-381 CWV-7dpi-1 GIM-361
## 4  GIM-382    g_Serratia     0 CqW-7dpi-2 GIM-382 CWV-7dpi-2 GIM-362
## 5  GIM-382    g_Bacillus     3 CqW-7dpi-2 GIM-382 CWV-7dpi-2 GIM-362
## 6  GIM-382 g_Pseudomonas   110 CqW-7dpi-2 GIM-382 CWV-7dpi-2 GIM-362
##         Index        mosq_species Treatment  dpi group.plaque.qPCR.1000.
## 1 CWV-7dpi-B1 Cx_quinquefasciatus    CQ-WNV 7dpi                  WNV(-)
## 2 CWV-7dpi-B1 Cx_quinquefasciatus    CQ-WNV 7dpi                  WNV(-)
## 3 CWV-7dpi-B1 Cx_quinquefasciatus    CQ-WNV 7dpi                  WNV(-)
## 4 CWV-7dpi-B2 Cx_quinquefasciatus    CQ-WNV 7dpi                  WNV(-)
## 5 CWV-7dpi-B2 Cx_quinquefasciatus    CQ-WNV 7dpi                  WNV(-)
## 6 CWV-7dpi-B2 Cx_quinquefasciatus    CQ-WNV 7dpi                  WNV(-)
##   group.plaque.qPCR.0. group.plaque.qPCR.1000.body.1000.
## 1                                                       
## 2                                                       
## 3                                                       
## 4                                                       
## 5                                                       
## 6                                                       
##   group.plaque.qPCR.0.body.0.       Info1 plaque_assay qPCR_head._ZIKV_WNV
## 1               7dpi-WNV(-/+) CQ-7dpi-WNV  7dpi-WNV(-)                   0
## 2               7dpi-WNV(-/+) CQ-7dpi-WNV  7dpi-WNV(-)                   0
## 3               7dpi-WNV(-/+) CQ-7dpi-WNV  7dpi-WNV(-)                   0
## 4               7dpi-WNV(-/+) CQ-7dpi-WNV  7dpi-WNV(-)                   0
## 5               7dpi-WNV(-/+) CQ-7dpi-WNV  7dpi-WNV(-)                   0
## 6               7dpi-WNV(-/+) CQ-7dpi-WNV  7dpi-WNV(-)                   0
##   qPCR_body._ZIKV_WNV head_ZIKV.1000 body_ZIKV.1000 Sample_or_Control    conc
## 1                2399                                            TRUE 21.9990
## 2                2399                                            TRUE 21.9990
## 3                2399                                            TRUE 21.9990
## 4               16184                                            TRUE 10.7685
## 5               16184                                            TRUE 10.7685
## 6               16184                                            TRUE 10.7685
df_merge$log10reads <- log10(df_merge$value)
is.na(df_merge$log10reads) <- sapply(df_merge$log10reads, is.infinite)
df_merge$log10reads[is.na(df_merge$log10reads)] <- 0
df_merge$Genus <- gsub("g_", "", df_merge$Genus)
df_merge$Genus <- factor(df_merge$Genus, levels = order); levels(df_merge$Genus)
## [1] "Serratia"    "Pseudomonas" "Bacillus"
df_merge_cq <- df_merge
df_merge_cq <- df_merge_cq[,colnames(df_merge_cq) %in% c("Treatment", "Genus", "variable", "log10reads")]

#### Aedes & Culex ####
df_merge_all <- rbind(df_merge_aa, df_merge_cq)
df_merge_all$Treatment
##    [1] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##    [7] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##   [13] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##   [19] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##   [25] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##   [31] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##   [37] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##   [43] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##   [49] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##   [55] AA-blood   AA-blood   AA-blood   AA-blood   AA-blood   AA-blood  
##   [61] AA-blood   AA-blood   AA-blood   AA-blood   AA-blood   AA-blood  
##   [67] AA-blood   AA-blood   AA-blood   AA-blood   AA-blood   AA-blood  
##   [73] AA-blood   AA-blood   AA-blood   AA-blood   AA-blood   AA-blood  
##   [79] AA-blood   AA-blood   AA-blood   AA-ZIKV    AA-ZIKV    AA-ZIKV   
##   [85] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##   [91] AA-blood   AA-blood   AA-blood   AA-blood   AA-blood   AA-blood  
##   [97] AA-blood   AA-blood   AA-blood   AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [103] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [109] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [115] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-blood   AA-blood   AA-blood  
##  [121] AA-blood   AA-blood   AA-blood   AA-blood   AA-blood   AA-blood  
##  [127] AA-blood   AA-blood   AA-blood   AA-blood   AA-blood   AA-blood  
##  [133] AA-blood   AA-blood   AA-blood   AA-blood   AA-blood   AA-blood  
##  [139] AA-blood   AA-blood   AA-blood   AA-blood   AA-blood   AA-blood  
##  [145] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [151] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-blood   AA-blood   AA-blood  
##  [157] AA-blood   AA-blood   AA-blood   AA-blood   AA-blood   AA-blood  
##  [163] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [169] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [175] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [181] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
##  [187] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
##  [193] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
##  [199] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
##  [205] AA-sucrose AA-sucrose AA-sucrose AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [211] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [217] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
##  [223] AA-sucrose AA-sucrose AA-sucrose AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [229] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [235] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [241] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-sucrose AA-sucrose AA-sucrose
##  [247] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
##  [253] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
##  [259] AA-sucrose AA-sucrose AA-sucrose AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [265] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [271] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
##  [277] AA-sucrose AA-sucrose AA-sucrose AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [283] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [289] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [295] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [301] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [307] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [313] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [319] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [325] AA-blood   AA-blood   AA-blood   AA-blood   AA-blood   AA-blood  
##  [331] AA-blood   AA-blood   AA-blood   AA-blood   AA-blood   AA-blood  
##  [337] AA-blood   AA-blood   AA-blood   AA-blood   AA-blood   AA-blood  
##  [343] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
##  [349] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
##  [355] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
##  [361] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
##  [367] AA-sucrose AA-sucrose AA-sucrose AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [373] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [379] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [385] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [391] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [397] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [403] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [409] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [415] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [421] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [427] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [433] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [439] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [445] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [451] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [457] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [463] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [469] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [475] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [481] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [487] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [493] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [499] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [505] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [511] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [517] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [523] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [529] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [535] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [541] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [547] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [553] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [559] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [565] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [571] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [577] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [583] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [589] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [595] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [601] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [607] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [613] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [619] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [625] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [631] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [637] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [643] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [649] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [655] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [661] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [667] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [673] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [679] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [685] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [691] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [697] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [703] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [709] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [715] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [721] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [727] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [733] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [739] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [745] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [751] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [757] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [763] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [769] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [775] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [781] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [787] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [793] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [799] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [805] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [811] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [817] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [823] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [829] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [835] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [841] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [847] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [853] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [859] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [865] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [871] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [877] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [883] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [889] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [895] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [901] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [907] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [913] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [919] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [925] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [931] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [937] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [943] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [949] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [955] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [961] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [967] CQ-blood   CQ-blood   CQ-blood   CQ-blood   CQ-blood   CQ-blood  
##  [973] CQ-blood   CQ-blood   CQ-blood   CQ-blood   CQ-blood   CQ-blood  
##  [979] CQ-blood   CQ-blood   CQ-blood   CQ-blood   CQ-blood   CQ-blood  
##  [985] CQ-blood   CQ-blood   CQ-blood   CQ-blood   CQ-blood   CQ-blood  
##  [991] CQ-blood   CQ-blood   CQ-blood   CQ-blood   CQ-blood   CQ-blood  
##  [997] CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose
## [1003] CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose
## [1009] CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose
## [1015] CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose
## [1021] CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose
## Levels: AA-sucrose AA-blood AA-ZIKV CQ-sucrose CQ-blood CQ-WNV
df_merge_all$Treatment <- factor(df_merge_all$Treatment , levels = c("AA-sucrose", "AA-blood", "AA-ZIKV", "CQ-sucrose", "CQ-blood", "CQ-WNV"))
####  plot  ####
ggplot(df_merge_all, aes(x = Treatment, y = log10reads, fill = Treatment)) + 
  facet_wrap(~Genus,ncol = 3, scales="free_x") + 
  theme_bw() +
  geom_boxplot(na.rm = TRUE,outlier.shape = NA, alpha = 0.7)+
  geom_jitter(shape=19, position=position_jitter(0.2), size = 2,alpha = 0.5)+
  scale_fill_manual(values = c(rev(col_6_Ae[c(4:6)]), rev(col_6_Ae[c(4:6)])))+
  theme(axis.title.y = element_text(size = 14, vjust=0.5),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        legend.position = "right",
        strip.text.x =element_text(size = 14, hjust = 0.5)) +
  ylab("log10 rarefied counts")

5 Phageome analysis

GuadPhage <- readRDS("~/Desktop/mosq_Guad_infection/RDS/GIM_phagome_raw.RDS")

5.1 HEATMAP of phageome

#### load phyloseq project #### 
Phage_host = readRDS("~/Desktop/mosq_Guad_infection/RDS/GIM_phage_host.RDS")
df <- as.data.frame(otu_table(Phage_host))
df_tax <- as.data.frame(tax_table(Phage_host))
df_meta <- data.frame(sample_data(Phage_host))

cat <- c("AA-sucrose", "AA-blood", "AA-ZIKV",  "CQ-sucrose", "CQ-blood", "CQ-WNV")
#### mat1 log10 reads per treatment #####
Phage_host_treatment <- Phage_host %>%
  tax_glom(taxrank = "Treatment")

Phage_host_treatment_merge <- merge(as.data.frame(Phage_host_treatment@otu_table), as.data.frame(Phage_host_treatment@tax_table), by="row.names")
rownames(Phage_host_treatment_merge) <- Phage_host_treatment_merge$Treatment
Phage_host_treatment_merge <- Phage_host_treatment_merge[, !colnames(Phage_host_treatment_merge) %in% c("Row.names", "mosq_species", "Treatment","dpi", "Info1")] 
mat1 <- log10(Phage_host_treatment_merge)
is.na(mat1) <- sapply(mat1, is.infinite)
mat1[is.na(mat1)]<-0 
mat1[1:5,1:5]
##            GIM319_NODE_68_length_3999_cov_6.501275
## AA-blood                                  3.611723
## AA-sucrose                                3.623559
## CQ-WNV                                    0.000000
## CQ-blood                                  0.000000
## CQ-sucrose                                0.000000
##            GIM424_NODE_178_length_2803_cov_5.668379
## AA-blood                                   2.209515
## AA-sucrose                                 0.000000
## CQ-WNV                                     3.046105
## CQ-blood                                   0.000000
## CQ-sucrose                                 0.000000
##            GIM573_NODE_150_length_3346_cov_8.985317
## AA-blood                                   3.247237
## AA-sucrose                                 3.092370
## CQ-WNV                                     0.000000
## CQ-blood                                   0.000000
## CQ-sucrose                                 0.000000
##            GIM95_NODE_22_length_3970_cov_5.511174
## AA-blood                                 3.411620
## AA-sucrose                               3.232742
## CQ-WNV                                   0.000000
## CQ-blood                                 0.000000
## CQ-sucrose                               0.000000
##            GIM400_NODE_1_length_17961_cov_34.782767
## AA-blood                                   2.460898
## AA-sucrose                                 0.000000
## CQ-WNV                                     6.087870
## CQ-blood                                   5.796710
## CQ-sucrose                                 4.686609
mat1 <- mat1[cat,]
rownames(mat1)
## [1] "AA-sucrose" "AA-blood"   "AA-ZIKV"    "CQ-sucrose" "CQ-blood"  
## [6] "CQ-WNV"
#### mat2 presence% per treatment #####
Phage_host_pre <- Phage_host %>%
  tax_glom(taxrank = "Treatment")
Phage_host_pre_merge <- merge(as.data.frame(Phage_host_pre@otu_table), as.data.frame(Phage_host_pre@tax_table), by="row.names")
rownames(Phage_host_pre_merge) <- Phage_host_pre_merge$Treatment
Phage_host_pre_merge <- Phage_host_pre_merge[, !colnames(Phage_host_pre_merge) %in% c("Row.names", "mosq_species", "Treatment","dpi", "Info1")] 
mat2 <- Phage_host_pre_merge
mat2 <- mat2[cat,]
rownames(mat2)
## [1] "AA-sucrose" "AA-blood"   "AA-ZIKV"    "CQ-sucrose" "CQ-blood"  
## [6] "CQ-WNV"
table(df_tax$Treatment)
## 
##   AA-blood AA-sucrose    AA-ZIKV   CQ-blood CQ-sucrose     CQ-WNV 
##         10         10         74         10         10         40
mat2[1,] <- mat2[1,]/10*100
mat2[2,] <- mat2[2,]/10*100
mat2[3,] <- mat2[3,]/74*100
mat2[4,] <- mat2[4,]/10*100
mat2[5,] <- mat2[5,]/10*100
mat2[6,] <- mat2[6,]/40*100


##abundance log transform
dflog <- log10(df)
is.na(dflog) <- sapply(dflog, is.infinite)
dflog[is.na(dflog)]<-0

#### order rows-sample and columns-contigs #### 
##column - phage host
host_or <- data.frame(sort(table(df_meta$host_genus), decreasing=T)); host_or
##           Var1 Freq
## 1    no result  203
## 2    Wolbachia   23
## 3 Enterobacter   21
## 4     Serratia    8
## 5  Phytobacter    6
## 6        Asaia    2
df_meta$host_genus <- factor(df_meta$host_genus, levels = host_or$Var1); levels(df_meta$host_genus)
## [1] "no result"    "Wolbachia"    "Enterobacter" "Serratia"     "Phytobacter" 
## [6] "Asaia"
##column - viral annotation
family_or <- data.frame(sort(table(df_meta$family), decreasing=T)); family_or
##             Var1 Freq
## 1  no annotation  222
## 2   Siphoviridae   17
## 3     Myoviridae    9
## 4    Podoviridae    7
## 5   Microviridae    6
## 6 Herelleviridae    2
df_meta$family <- factor(df_meta$family, levels = family_or$Var1); levels(df_meta$family)
## [1] "no annotation"  "Siphoviridae"   "Myoviridae"     "Podoviridae"   
## [5] "Microviridae"   "Herelleviridae"
##column - sample order-length
df_meta$log10length <- log10(df_meta$length)
contigs_ordered <- df_meta[with(df_meta, order(host_genus, -log10length)),]
df_meta$length1k <- df_meta$length/1000

########  draw the map ########
heatmap_col_tp1 <- colorRamp2(c(0:5), brewer.pal(6, "YlOrBr"), space = "RGB") #heatmap colour
heatmap_col_tp2 <- colorRamp2(c(0,20,40,60,80,100), brewer.pal(6, "YlGnBu"), space = "RGB")
heatmap_col_tp3 <- colorRamp2(c(0,20,40,60,80,100), brewer.pal(6, "Greens"), space = "RGB")

#### host block ####  
colh = c("#818687FF", "#2D5C21FF", "#2CA0A9FF", "#A1E4BDFF", "#323B34FF", "#C8D656FF")
colhost = HeatmapAnnotation(host = anno_block(gp = gpar(fill = colh)),show_annotation_name = F, annotation_label = F) 

#### viral annotation ####  
colf = c("#818687FF", "#66C2A5", "#FC8D62", "#8DA0CB", "#FFD92F", "#E5C494")
colfamily = HeatmapAnnotation(condition = df_meta$family,
                              col = list(condition = c("no annotation"=colf[1], "Siphoviridae"=colf[2], "Myoviridae"=colf[3],
                                                       "Podoviridae"=colf[4], "Microviridae"=colf[5], "Herelleviridae"=colf[6]),
                                                       show_annotation_name = T, annotation_label = F, 
                                                       border = T))

#### contigs length ####  
collen = HeatmapAnnotation(length = anno_barplot(df_meta$length1k, ylim = c(0,50),
                                                 height = unit(2, "cm")))
df_meta$aai_completeness[is.na(df_meta$aai_completeness)] <- 0
complete = HeatmapAnnotation(complete = df_meta$aai_completeness, col = list(complete=heatmap_col_tp3))

#### legend ####  
lgd_host = Legend(labels = levels(df_meta$host_genus), title = "Predicted_host", legend_gp = gpar(fill = colh))
lgd_virus = Legend(labels = levels(df_meta$family), title = "Phage_annotation", legend_gp = gpar(fill = colf))
# lgd_reads = Legend(title = "Log10 reads number", col = heatmap_col_tp1, at = c(0:6), labels = c("0", "1", "2", "3", "4", "5", "6"))
# lgd_pro = Legend(title = "Relative presence", col = heatmap_col_tp2, at = c(0:5), labels = c("0", "20", "40", "60", "80", "100"))
pd = packLegend(lgd_host, lgd_virus, 
                # lgd_reads, lgd_pro, 
                direction = "vertical")
#### plot ####  
ht1 <- Heatmap(as.matrix(mat1),
               name = "log10 reads number", #title of legend
               row_title = "Samples",
               col = heatmap_col_tp1,
               column_split = df_meta$host_genus, 
               row_names_gp = gpar(fontsize = 12),
               column_names_gp = gpar(fontsize = 0),
               column_order = contigs_ordered$name,
               row_order = cat,
               cluster_rows = F,
               cluster_columns= F,
               rect_gp = gpar(col = "white", lwd = 0.5))
ht2 <- Heatmap(as.matrix(mat2),
               name = "Present proportion", #title of legend
               row_title = "Samples",
               col = heatmap_col_tp2,
               column_split = df_meta$host_genus, 
               row_names_gp = gpar(fontsize = 12),
               column_names_gp = gpar(fontsize = 0),
               column_order = contigs_ordered$name,
               row_order = cat,
               cluster_rows = F,
               cluster_columns= F,
               rect_gp = gpar(col = "white", lwd = 0.5))
 
ht_list = colhost %v% colfamily  %v% complete %v% collen %v% ht1 %v% ht2
# draw(ht_list, annotation_legend_list = pd)
ht_draw = draw(ht_list, annotation_legend_list = pd)

5.2 Alpha diversity of phageome

GuadPhage_samselect = prune_samples(sample_sums(GuadPhage) > 0 , GuadPhage)
GuadPhage_samselect
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 263 taxa and 150 samples ]
## sample_data() Sample Data:       [ 150 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 263 taxa by 17 taxonomic ranks ]
factor(sample_data(GuadPhage_samselect)$Treatment)
##   [1] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##   [7] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [13] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [19] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [25] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [31] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [37] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [43] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [49] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [55] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [61] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV   
##  [67] AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-ZIKV    AA-sucrose
##  [73] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
##  [79] AA-sucrose AA-sucrose AA-sucrose AA-blood   AA-blood   AA-blood  
##  [85] AA-blood   AA-blood   AA-blood   AA-blood   AA-blood   AA-blood  
##  [91] AA-blood   CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
##  [97] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
## [103] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
## [109] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
## [115] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
## [121] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV    
## [127] CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-WNV     CQ-sucrose
## [133] CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose
## [139] CQ-sucrose CQ-sucrose CQ-blood   CQ-blood   CQ-blood   CQ-blood  
## [145] CQ-blood   CQ-blood   CQ-blood   CQ-blood   CQ-blood   CQ-blood  
## Levels: AA-blood AA-sucrose AA-ZIKV CQ-blood CQ-sucrose CQ-WNV
sample_data(GuadPhage_samselect)$Treatment <- factor(sample_data(GuadPhage_samselect)$Treatment, 
                                                     levels = c("AA-sucrose","AA-blood", "AA-ZIKV", "CQ-sucrose", "CQ-blood", "CQ-WNV"))

phage <- scan(file = "~/OneDrive/1. Project + Guadeloupe mosq/phage_contigs_complete+20.txt", what="", sep="\n")
my_subset <- subset(otu_table(GuadPhage_samselect), rownames(otu_table(GuadPhage_samselect)) %in% phage)
new_GuadPhage_samselect <- merge_phyloseq(my_subset, tax_table(GuadPhage_samselect), sample_data(GuadPhage_samselect))

plot_richness(new_GuadPhage_samselect, x="Treatment",measures=c("Shannon"),col="Treatment", shape = "mosq_species") + 
  scale_color_manual(values = rev(c(col[9],col[5],col[2],col[9],col[5],col[2]))) +
  theme_bw()+
  # facet_wrap(~mosq_species, scales = "free")+
  scale_shape_manual(values = c(16, 17))+
  theme(axis.title.x = element_blank(),
        axis.text.x  = element_blank(),
        strip.text.x = element_text(size = 12, colour = "black"),
        plot.title = element_text(hjust = 0.5, size = 12),
        axis.ticks.x = element_blank()) +
  ggtitle("Alpha diversity of phages")+
  ylab("Alpha Diversity Value")
## Warning in estimate_richness(physeq, split = TRUE, measures = measures): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
## 
## We recommended that you find the un-trimmed data and retry.

erich <- estimate_richness(new_GuadPhage_samselect, measures = c("Shannon"))
## Warning in estimate_richness(new_GuadPhage_samselect, measures = c("Shannon")): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
## 
## We recommended that you find the un-trimmed data and retry.
pairwise.wilcox.test(erich$Shannon,sample_data(new_GuadPhage_samselect)$Treatment, p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  erich$Shannon and sample_data(new_GuadPhage_samselect)$Treatment 
## 
##            AA-sucrose AA-blood AA-ZIKV CQ-sucrose CQ-blood
## AA-blood   1.00000    -        -       -          -       
## AA-ZIKV    1.00000    1.00000  -       -          -       
## CQ-sucrose 0.00051    0.00051  1.1e-08 -          -       
## CQ-blood   0.00026    0.00022  5.3e-10 0.08898    -       
## CQ-WNV     4.1e-06    4.1e-06  < 2e-16 0.00047    0.14860 
## 
## P value adjustment method: BH

5.3 Phageome of Aedes aegypti

####  color ####
col = plasma(9,alpha = 1, begin=0.9, end = 0, direction = 1)

#### Phyloseq Project #### 
GuadPhage_Ae_aegypti <- subset_samples(GuadPhage, mosq_species  == "Ae_aegypti" ) 
GuadPhage_Ae_aegypti
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 263 taxa and 94 samples ]
## sample_data() Sample Data:       [ 94 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 263 taxa by 17 taxonomic ranks ]
GP.M1_AA = GuadPhage_Ae_aegypti
GP.M1_AA
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 263 taxa and 94 samples ]
## sample_data() Sample Data:       [ 94 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 263 taxa by 17 taxonomic ranks ]
sum(sample_sums(GP.M1_AA))
## [1] 182924
sample_data(GP.M1_AA)$Treatment <- factor(sample_data(GP.M1_AA)$Treatment, levels = c("AA-ZIKV","AA-blood","AA-sucrose"))

#### subset #### 
sample_sums(GP.M1_AA)
##  GIM21  GIM22  GIM23  GIM24  GIM25  GIM26  GIM27  GIM28  GIM29  GIM30  GIM31 
##   1212   1197   1284    359   1381   1432   1242   1345   1070   1418   1386 
##  GIM32  GIM33  GIM34  GIM35  GIM36  GIM37  GIM38  GIM39  GIM40  GIM81  GIM82 
##   1633   1047    602   1881   1554   1280   1405   1067   1382   1899    996 
##  GIM83  GIM84  GIM85  GIM86  GIM87  GIM88  GIM89  GIM90  GIM91  GIM92  GIM93 
##    330    467   1048    719   1043    438    175   1841    927      0    981 
##  GIM94  GIM95  GIM96  GIM97  GIM98  GIM99 GIM100 GIM572 GIM573 GIM574 GIM575 
##    897    823   1723   1425   1426    817   1418   1317   2604   1496    971 
## GIM576 GIM577 GIM578 GIM580 GIM581 GIM582 GIM583 GIM585 GIM586 GIM608 GIM609 
##   2124   1525   2095    854   1038      0   2319   1261   1389   1470    598 
## GIM610 GIM611 GIM612 GIM613 GIM614 GIM615 GIM616 GIM617 GIM618 GIM619 GIM620 
##   1299   1229    434   1062   1176   1251    525      0   2383   1462    586 
## GIM621 GIM622 GIM623 GIM624 GIM625 GIM626 GIM627 GIM628 GIM317 GIM318 GIM319 
##    848    919    909    666    841   1550   1432    410    732    539   1637 
## GIM320 GIM347 GIM349 GIM350 GIM549 GIM553 GIM554 GIM257 GIM258 GIM259 GIM260 
##   1688   1339   1088  77045   2418   1118   1259    654   2376    722    456 
## GIM284 GIM287 GIM288 GIM290 GIM543 GIM546 
##   1188   1686    483    923    150    810
sample_sums(GP.M1_AA)[which(sample_sums(GP.M1_AA) == 0)]
##  GIM92 GIM582 GIM617 
##      0      0      0
GP.M1_AA_samselect = prune_samples(sample_sums(GP.M1_AA) > 0 , GP.M1_AA)
GP.M1_AA_samselect
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 263 taxa and 91 samples ]
## sample_data() Sample Data:       [ 91 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 263 taxa by 17 taxonomic ranks ]
GP.M1_AA_select_0 = prune_taxa(taxa_sums(GP.M1_AA_samselect) > 0, GP.M1_AA_samselect)
GP.M1_AA_select_0
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 151 taxa and 91 samples ]
## sample_data() Sample Data:       [ 91 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 151 taxa by 17 taxonomic ranks ]
#### PCoA ####
GP.ord.GP.M1_AA_samselect <- ordinate(GP.M1_AA_samselect, "PCoA", "bray")
plot_ordination(GP.M1_AA_samselect, GP.ord.GP.M1_AA_samselect, type = "samples", col = "Treatment")+#,label= "name") +
  theme_bw() +
  geom_point(size = 4, alpha=0.7)  + 
  scale_color_manual(values = c(col[9],col[5],col[2]))+
  theme(legend.title=element_blank(),
        legend.position = "right",
        plot.title = element_text(size = 12,hjust = 0.5))+
  ggtitle("PCoA of phages in all Ae. aegypti samples")

set.seed(1)
# Calculate bray curtis distance matrix
GP.M1_species_select_bray <- phyloseq::distance(GP.M1_AA_samselect, method = "bray")
# make a data frame from the sample_data
GP.M1_species_select_sampledf <- data.frame(sample_data(GP.M1_AA_samselect))
# Adonis test
adonis(GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf)
## 
## Call:
## adonis(formula = GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)  
## Treatment  2    0.7636 0.38180  2.2902 0.04947  0.027 *
## Residuals 88   14.6707 0.16671         0.95053         
## Total     90   15.4343                 1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5.4 Phageome Culex quinquefasciatus

#### Phyloseq Project #### 
GuadPhage_Cx_quinquefasciatus <- subset_samples(GuadPhage, mosq_species  == "Cx_quinquefasciatus" ) 
GuadPhage_Cx_quinquefasciatus
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 263 taxa and 60 samples ]
## sample_data() Sample Data:       [ 60 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 263 taxa by 17 taxonomic ranks ]
GP.M1_CQ = GuadPhage_Cx_quinquefasciatus
GP.M1_CQ
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 263 taxa and 60 samples ]
## sample_data() Sample Data:       [ 60 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 263 taxa by 17 taxonomic ranks ]
sample_data(GP.M1_CQ)$group.plaque.qPCR.0.body.0. <- factor(sample_data(GP.M1_CQ)$group.plaque.qPCR.0.body.0., 
                                                            levels = c("7dpi-WNV(+/)", "7dpi-WNV(-/+)", "7dpi-WNV(-/-)","CQ-7dpi-blood","CQ-7dpi-sucrose",
                                                                       "14dpi-WNV(+/)","14dpi-WNV(-/+)","14dpi-WNV(-/-)","CQ-14dpi-blood","CQ-14dpi-sucrose"))
sample_data(GP.M1_CQ)$Info1 <- factor(sample_data(GP.M1_CQ)$Info1, levels = c("CQ-7dpi-WNV", "CQ-7dpi-blood", "CQ-7dpi-sucrose", "CQ-14dpi-WNV", "CQ-14dpi-blood", "CQ-14dpi-sucrose"))
sample_data(GP.M1_CQ)$Treatment <- factor(sample_data(GP.M1_CQ)$Treatment, levels = c("CQ-WNV","CQ-blood","CQ-sucrose"))
levels(sample_data(GP.M1_CQ)$Info1)
## [1] "CQ-7dpi-WNV"      "CQ-7dpi-blood"    "CQ-7dpi-sucrose"  "CQ-14dpi-WNV"    
## [5] "CQ-14dpi-blood"   "CQ-14dpi-sucrose"
levels(sample_data(GP.M1_CQ)$group.plaque.qPCR.0.body.0.)
##  [1] "7dpi-WNV(+/)"     "7dpi-WNV(-/+)"    "7dpi-WNV(-/-)"    "CQ-7dpi-blood"   
##  [5] "CQ-7dpi-sucrose"  "14dpi-WNV(+/)"    "14dpi-WNV(-/+)"   "14dpi-WNV(-/-)"  
##  [9] "CQ-14dpi-blood"   "CQ-14dpi-sucrose"
#### subset #### 
sample_sums(GP.M1_CQ)
## GIM381 GIM382 GIM383 GIM384 GIM385 GIM386 GIM387 GIM388 GIM389 GIM390 GIM391 
##  32527  48435 106963 180814 115000 223177 223468 197054 108111 420751  66868 
## GIM392 GIM393 GIM394 GIM395 GIM396 GIM397 GIM398 GIM399 GIM400 GIM421 GIM422 
##  14533 717060 216103 262323 389183 129958 153894 423605  35159 406227  80616 
## GIM423 GIM424 GIM425 GIM426 GIM427 GIM428 GIM429 GIM430 GIM431 GIM432 GIM433 
## 363912 114523 135616  83275  71798 177505  59004  83063  97328  84743 130900 
## GIM434 GIM435 GIM436 GIM437 GIM438 GIM439 GIM440 GIM491 GIM492 GIM493 GIM494 
##  43898 138201  36543  91653  98097 736556  17718 153484  35892   2698   6944 
## GIM495 GIM525 GIM527 GIM528 GIM529 GIM530 GIM451 GIM452 GIM453 GIM454 GIM455 
##   4822  26017  10736  79978  21879      0 574667 194018  87396 536174  71076 
## GIM471 GIM472 GIM473 GIM475 GIM476 
## 479370 572072   7832 369564 599136
sample_sums(GP.M1_CQ)[which(sample_sums(GP.M1_CQ) == 0)]
## GIM530 
##      0
GP.M1_CQ_samselect = prune_samples(sample_sums(GP.M1_CQ) > 0 , GP.M1_CQ)
GP.M1_CQ_samselect
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 263 taxa and 59 samples ]
## sample_data() Sample Data:       [ 59 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 263 taxa by 17 taxonomic ranks ]
GP.M1_CQ_select_0 = prune_taxa(taxa_sums(GP.M1_CQ_samselect) > 0, GP.M1_CQ_samselect)
GP.M1_CQ_select_0
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 164 taxa and 59 samples ]
## sample_data() Sample Data:       [ 59 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 164 taxa by 17 taxonomic ranks ]
#### PCoA all samples #### 
GP.ord.GP.M1_CQ_samselect <- ordinate(GP.M1_CQ_samselect, "PCoA", "bray")

plot_ordination(GP.M1_CQ_samselect, GP.ord.GP.M1_CQ_samselect, type = "samples", col = "Treatment")+
  theme_bw() +
  geom_point(size = 4, alpha=0.7)  + 
  scale_color_manual(values = c(col[9],col[5],col[2]))+
  theme(legend.title=element_blank(),
        legend.position = "right",
        plot.title = element_text(size = 12,hjust = 0.5))+
  ggtitle("PCoA of phages in all Cx. quinquefasciatus samples")+
  stat_ellipse(type = "norm", linetype = 2,level=0.6,aes(group = Treatment))

set.seed(1)
# Calculate bray curtis distance matrix
GP.M1_species_select_bray <- phyloseq::distance(GP.M1_CQ_select_0, method = "bray")
# make a data frame from the sample_data
GP.M1_species_select_sampledf <- data.frame(sample_data(GP.M1_CQ_select_0))
# Adonis test
adonis(GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf)
## 
## Call:
## adonis(formula = GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Treatment  2    2.5667 1.28336  6.0101 0.17671  0.001 ***
## Residuals 56   11.9580 0.21353         0.82329           
## Total     58   14.5247                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# # Call:
beta <- betadisper(GP.M1_species_select_bray, GP.M1_species_select_sampledf$Treatment)
permutest(beta)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df Sum Sq  Mean Sq      F N.Perm Pr(>F)
## Groups     2 0.0167 0.008348 0.2089    999  0.783
## Residuals 56 2.2379 0.039962
#### PCoA 7 dpe #### 
GP.M1_CQ_select_0_7dpi = subset_samples(GP.M1_CQ_samselect, dpi == "7dpi")
GP.M1_CQ_select_7dpi = prune_taxa(taxa_sums(GP.M1_CQ_select_0_7dpi) > 0, GP.M1_CQ_select_0_7dpi)
levels(sample_data(GP.M1_CQ_select_7dpi)$group.plaque.qPCR.0.body.0.)
## [1] "7dpi-WNV(+/)"    "7dpi-WNV(-/+)"   "7dpi-WNV(-/-)"   "CQ-7dpi-blood"  
## [5] "CQ-7dpi-sucrose"
GP.ord.GP.M1_CQ_select_7dpi <- ordinate(GP.M1_CQ_select_7dpi, "PCoA", "bray")

plot_ordination(GP.M1_CQ_select_7dpi, GP.ord.GP.M1_CQ_select_7dpi, type = "samples", color = "Treatment", shape = "group.plaque.qPCR.0.body.0.") +
  theme_bw()+
  geom_point(size = 4, alpha=0.7)  + 
  scale_shape_manual(values=c(shap[1:4],shap[6])) +
  scale_color_manual(values = c(col[9],col[5],col[2]))+
  theme(legend.title=element_blank(),
        legend.position = "right",
        plot.title = element_text(size = 12,hjust = 0.5))+
  ggtitle("PCoA of phages in Cx. quinquefasciatus at 7 dpe")+
  stat_ellipse(type = "norm", linetype = 2,level=0.6,aes(group = Treatment), color="black")

set.seed(1)
# Calculate bray curtis distance matrix
GP.M1_species_select_bray <- phyloseq::distance(GP.M1_CQ_select_0_7dpi, method = "bray")
# make a data frame from the sample_data
GP.M1_species_select_sampledf <- data.frame(sample_data(GP.M1_CQ_select_0_7dpi))
# Adonis test
adonis(GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf)
## 
## Call:
## adonis(formula = GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## Treatment  2    1.2973 0.64865  2.8944 0.17655  0.005 **
## Residuals 27    6.0509 0.22411         0.82345          
## Total     29    7.3482                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#### PCoA 14 dpe #### 
GP.M1_CQ_select_0_14dpi = subset_samples(GP.M1_CQ_samselect, dpi == "14dpi")
GP.M1_CQ_select_14dpi = prune_taxa(taxa_sums(GP.M1_CQ_select_0_14dpi) > 0, GP.M1_CQ_select_0_14dpi)
levels(sample_data(GP.M1_CQ_select_14dpi)$group.plaque.qPCR.0.body.0.)
## [1] "14dpi-WNV(+/)"    "14dpi-WNV(-/+)"   "14dpi-WNV(-/-)"   "CQ-14dpi-blood"  
## [5] "CQ-14dpi-sucrose"
GP.ord.GP.M1_CQ_select_14dpi <- ordinate(GP.M1_CQ_select_14dpi, "PCoA", "bray")

plot_ordination(GP.M1_CQ_select_14dpi, GP.ord.GP.M1_CQ_select_14dpi, type = "samples", color = "Treatment", shape = "group.plaque.qPCR.0.body.0.")+
  theme_bw()+
  geom_point(size = 4, alpha=0.7)  +  
  scale_shape_manual(values=c(shap[1:4],shap[6])) +
  scale_color_manual(values = c(col[9],col[5],col[2]))+
  theme(legend.title=element_blank(),
        legend.position = "right",
        plot.title = element_text(size = 12,hjust = 0.5))+
  ggtitle("PCoA of phages in Cx. quinquefasciatus at 14dpi")+
  stat_ellipse(type = "norm", linetype = 2,level=0.8,aes(group = Treatment), color="black")

set.seed(1)
# Calculate bray curtis distance matrix
GP.M1_species_select_bray <- phyloseq::distance(GP.M1_CQ_select_0_14dpi, method = "bray")
# make a data frame from the sample_data
GP.M1_species_select_sampledf <- data.frame(sample_data(GP.M1_CQ_select_0_14dpi))
# Adonis test
adonis(GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf)
## 
## Call:
## adonis(formula = GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Treatment  2    1.9906 0.99531  5.2418 0.28735  0.001 ***
## Residuals 26    4.9368 0.18988         0.71265           
## Total     28    6.9275                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5.5 Phage correlate with bacteria

####  phage  ####
Phage_host = readRDS("~/Desktop/mosq_Guad_infection/RDS/GIM_phage_host.RDS")
Phage_host_subset <- subset_samples(Phage_host, host_genus != "no result")
pOTU <- as.data.frame(t(Phage_host_subset@otu_table))
pTAX <- data.frame(sample_data(Phage_host_subset))
pTAX <- pTAX[,15:16]
pTAX$phageID <- rownames(pTAX)
pOTU$phageID <- rownames(pOTU)
phage_merge <- merge(pOTU, pTAX, by = "phageID")

####  bacteria  ####
Bac = readRDS("~/Desktop/mosq_Guad_infection/16s_rRNA/GIM_Bacteriome_decont_SLV132.rds")
BacGenus <- Bac %>%
  tax_glom(taxrank = "Genus")

BacOTU <- as.data.frame(t(BacGenus@otu_table))
colnames(BacOTU) <- paste("Bac_", colnames(BacOTU), sep = "")
BacTAX <- as.data.frame(BacGenus@tax_table)
BacOTU$bacID <- rownames(BacOTU)
BacTAX$bacID <- rownames(BacTAX)
Bac_merge <- merge(BacOTU, BacTAX[,5:7], by = "bacID")
Bac_merge$Genus <- gsub("g_", "", Bac_merge$Genus)

####  merge phage + bacteria  ####
colnames(phage_merge) <- gsub("host_genus", "Genus", colnames(phage_merge))

Bac_merge_slec <- Bac_merge[Bac_merge$Genus %in% c(phage_merge$Genus, "uc_f_Enterobacteriaceae"),]
phage_merge$Genus_og <- phage_merge$Genus
phage_merge$Genus <- gsub("Enterobacter", "uc_f_Enterobacteriaceae", phage_merge$Genus)
phage_merge$Genus <- gsub("Phytobacter", "uc_f_Enterobacteriaceae", phage_merge$Genus)

phage_Bac <- merge(phage_merge, Bac_merge_slec, all = T, by ="Genus") 
phage_Bac <- phage_Bac[, !colnames(phage_Bac) %in% c("Family", "host_species", "bacID", "Bac_NC-p07", "Bac_NC-p06", "Bac_GIM-NC1", "Bac_GIM-NC2")]
phage_Bac_melt <- melt(phage_Bac)
## Using Genus, phageID, Genus_og as id variables
phage_Bac_melt <- phage_Bac_melt %>% separate(variable, c("cat","Sample"), sep ="GIM"); head(phage_Bac_melt)
##      Genus                                   phageID Genus_og cat Sample value
## 1    Asaia GIM439_NODE_511_length_3345_cov_12.118421    Asaia         21     0
## 2    Asaia  GIM439_NODE_40_length_9507_cov_18.950795    Asaia         21     0
## 3 Serratia  GIM455_NODE_1_length_46241_cov_41.282255 Serratia         21     0
## 4 Serratia GIM476_NODE_105_length_6055_cov_20.624791 Serratia         21     0
## 5 Serratia GIM476_NODE_133_length_5420_cov_30.978664 Serratia         21     0
## 6 Serratia GIM476_NODE_406_length_3177_cov_20.609677 Serratia         21     0
phage_Bac_melt$Sample <- gsub("-", "", phage_Bac_melt$Sample)
phage_Bac_melt$Sample <- paste("GIM", phage_Bac_melt$Sample, sep = "")
phage_Bac_melt$cat <- ifelse(phage_Bac_melt$cat == "", "phage", "bacteria"); head(phage_Bac_melt)
##      Genus                                   phageID Genus_og   cat Sample
## 1    Asaia GIM439_NODE_511_length_3345_cov_12.118421    Asaia phage  GIM21
## 2    Asaia  GIM439_NODE_40_length_9507_cov_18.950795    Asaia phage  GIM21
## 3 Serratia  GIM455_NODE_1_length_46241_cov_41.282255 Serratia phage  GIM21
## 4 Serratia GIM476_NODE_105_length_6055_cov_20.624791 Serratia phage  GIM21
## 5 Serratia GIM476_NODE_133_length_5420_cov_30.978664 Serratia phage  GIM21
## 6 Serratia GIM476_NODE_406_length_3177_cov_20.609677 Serratia phage  GIM21
##   value
## 1     0
## 2     0
## 3     0
## 4     0
## 5     0
## 6     0
phage_Bac_melt_wide <- dcast(phage_Bac_melt, phageID + Sample + Genus + Genus_og ~ cat, value.var = 'value')

####  add meta  ####
head(phage_Bac_melt)
##      Genus                                   phageID Genus_og   cat Sample
## 1    Asaia GIM439_NODE_511_length_3345_cov_12.118421    Asaia phage  GIM21
## 2    Asaia  GIM439_NODE_40_length_9507_cov_18.950795    Asaia phage  GIM21
## 3 Serratia  GIM455_NODE_1_length_46241_cov_41.282255 Serratia phage  GIM21
## 4 Serratia GIM476_NODE_105_length_6055_cov_20.624791 Serratia phage  GIM21
## 5 Serratia GIM476_NODE_133_length_5420_cov_30.978664 Serratia phage  GIM21
## 6 Serratia GIM476_NODE_406_length_3177_cov_20.609677 Serratia phage  GIM21
##   value
## 1     0
## 2     0
## 3     0
## 4     0
## 5     0
## 6     0
phage_Bac_melt_wide$Genus <- as.factor(phage_Bac_melt_wide$Genus)

phage_Bac_melt_wide$log10bacteria <- log10(phage_Bac_melt_wide$bacteria)
is.na(phage_Bac_melt_wide$log10bacteria) <- sapply(phage_Bac_melt_wide$log10bacteria, is.infinite)
phage_Bac_melt_wide$log10bacteria[is.na(phage_Bac_melt_wide$log10bacteria)]<-0

phage_Bac_melt_wide$log10phage <- log10(phage_Bac_melt_wide$phage)
is.na(phage_Bac_melt_wide$log10phage) <- sapply(phage_Bac_melt_wide$log10phage, is.infinite)
phage_Bac_melt_wide$log10phage[is.na(phage_Bac_melt_wide$log10phage)]<-0

samplemeta <- data.frame(BacGenus@sam_data); colnames(samplemeta)
##  [1] "mosq"                              "body"                             
##  [3] "mosq_or"                           "head"                             
##  [5] "Index"                             "mosq_species"                     
##  [7] "Treatment"                         "dpi"                              
##  [9] "group.plaque.qPCR.1000."           "group.plaque.qPCR.0."             
## [11] "group.plaque.qPCR.1000.body.1000." "group.plaque.qPCR.0.body.0."      
## [13] "Info1"                             "plaque_assay"                     
## [15] "qPCR_head._ZIKV_WNV"               "qPCR_body._ZIKV_WNV"              
## [17] "head_ZIKV.1000"                    "body_ZIKV.1000"                   
## [19] "Sample_or_Control"                 "conc"
samplemeta <- samplemeta[ ,colnames(samplemeta) %in% c("mosq", "Index", "mosq_species", "Treatment", "dpi", "plaque_assay", "Info1")]
samplemeta$Sample <- gsub("-", "", rownames(samplemeta))

phage_Bac_melt_wide_meta <- merge(phage_Bac_melt_wide, samplemeta, by= "Sample")

levels(phage_Bac_melt_wide_meta$Genus)
## [1] "Asaia"                   "Serratia"               
## [3] "uc_f_Enterobacteriaceae" "Wolbachia"
phage_Bac_melt_wide_metaSle <- phage_Bac_melt_wide_meta[phage_Bac_melt_wide_meta$Genus %in% c("Wolbachia", "Serratia"),]

phage_Bac_melt_wide_metaSle$Genus <- factor(phage_Bac_melt_wide_metaSle$Genus, levels = c("Wolbachia", "Serratia"))
phage_Bac_melt_wide_metaSle$mosq_species <- factor(phage_Bac_melt_wide_metaSle$mosq_species, levels = c("Ae_aegypti", "Cx_quinquefasciatus"))
####  plot  ####
ggscatter(phage_Bac_melt_wide_metaSle, x = "log10bacteria", y = "log10phage",
          color = "mosq_species", palette = c("#7294D4", "#C27D38"), size = 2.2, alpha = 0.8,
          conf.int = T, cor.coef = TRUE, cor.method = "pearson")+
  facet_wrap(~Genus)+
  geom_smooth(color="black")+
  theme_bw()+
  theme(axis.title.y = element_text(size = 12, vjust=0.5),
        axis.title.x = element_text(size = 12, hjust=0.5),
        panel.grid.minor = element_blank(),
        legend.position = "top",
        strip.text.x =element_text(size = 14, hjust = 0.5)) +
  xlab("log10 host bacteria reads")+
  ylab("1og10 phage reads")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'